%  CalcSigVol.m
%     calculate thickness and volume of density layers
%
%          grid size in meters
dx=20d3;
dy=20d3;
xrange=[20 500];
yrange=[-20 300];
PlotNeg=1;

CruiseName='NBP0103';
%CruiseName='NBP0104';

load wap3_coast_xy.mat;

load([CruiseName '.gridstns.mat']);
load([CruiseName '.siggrid.mat']);

fid=fopen([CruiseName '.vols.prt'],'w');
fprintf(fid,[CruiseName '\n']);

fprintf(fid,'  sigma    sigma   volume (km^3)  avg thick (m)  npoints nnegthick\n');

[m nsig]=size(siglevels);
[n,nx]=size(xi);
[n,ny]=size(yi);
[xxi,yyi]=meshgrid(xi,yi);

disp(CruiseName)
disp('  sigma    sigma   volume (km^3)  avg thick (m)  npoints nnegthick')
for isig=1:nsig-1
  ThickSig=DepSigGrid(:,:,isig+1)-DepSigGrid(:,:,isig);
  ineg=find(ThickSig < 0);
  [negp n]=size(ineg);
%
  v=0;
  i=find(~isnan(ThickSig) );
  v=sum(ThickSig(i));
  [nv n]=size(i);
  vkm=1d-9*dx*dy*v;
  if(nv > 0) 
    avgth=v/nv;
  else
    avgth=0;
  end
  nvkm=0;
  if negp > 0
    v=sum(ThickSig(ineg));
     nvkm=1d-9*dx*dy*v;
  end

  disp([num2str(siglevels(isig+1)) ' ' num2str(siglevels(isig))  ...
        ' ' num2str(vkm) '  '  num2str(avgth) '  ' num2str(nv) '  ' ...
           num2str(negp) '  '  num2str(nvkm)]);
  disp('  ')
  fprintf(fid,'%8.3f %8.3f %11.4f %8.3f %8d %8d %11.4f\n',...
    siglevels(isig+1),siglevels(isig),vkm,avgth,nv,negp,nvkm);

  if (PlotNeg & negp > 0)
    figure
    plot(xcoast,ycoast);
    axis([xrange yrange]);
    xlabel('distance (km)');
    ylabel('distance (km)');
    title(['Points with negative thickness ' ...
           num2str(siglevels(isig+1)) ' ' num2str(siglevels(isig))]);
    hold on
    plot(xxi(ineg),yyi(ineg),'r*')
    hold off
    disp('  Negative thicknesses ')
     for j=1:negp
        disp([num2str(xxi(ineg(j))) ' ' num2str(yyi(ineg(j)))...
           '  ' num2str(ThickSig(ineg(j))) ])
     end
  end
end

fclose(fid);

