% CTDVIEW.M
% Interactive display of EPIC style netCDF CTD profiles.  
% Plots the station locations for the
% cruise, then lets the user click on the station location to display the 
% profiles at that station
%

% 2. Requires metcdf.mex, mcvgt.m, laxis.m, clearg.m, stations.m
%
% Rich Signell               |  rsignell@crusty.er.usgs.gov
% U.S. Geological Survey     |  (508) 457-2229  |  FAX (508) 457-2310
% Quissett Campus            | 
% Woods Hole, MA  02543      |
%--------------------------------------------------------------------

% the coastline file is hard-wired (mat file containing the variable 'coast')
clear
load mbay_coast

  var1='temp';
  var2='salt';
  var3='dens';
 
% stations.m  is an function that you need to create that returns a
% character matrix that contains the names of all the EPIC CTD profile
% files, one per line. 
% 
% EXAMPLE: function sta=stations();
%          sta=['mb06c001.cdf'
%               'mb06c002.cdf'
%               'mb06c003.cdf'];
%
%on UNIX: command=['! sed -f $HOME/m/ctdscript.m mb01.dat > tmpctd.m'];
%will create this file
%

sta=stations;

[m,n]=size(sta);
mexcdf('setopts',0);
for j=1:m
  cdfid=mexcdf('open',sta(j,:),'nowrite');
  lat(j)=mexcdf('varget1',cdfid,'latitude',0);
  lon(j)=-(mexcdf('varget1',cdfid,'longitude',0));
  rcode=mexcdf('close',cdfid);
end
%%plot the station locations
clearg
subplot(224)
axis('square')
ax=laxis(lon,lat);
mb=0;
axis(ax);
plot(lon,lat,'gx',coast(:,1),coast(:,2),'r-')
axis;
first=1;
i=sqrt(-1);
iinit=10;
icopy=iinit;
while(mb~=3),
  %%pick a station to plot
  text(.6,.03,'LEFT mouse on station','sc')
  text(.6,0,'RIGHT mouse to quit','sc')
  [x,y,mb]=ginput(1);
  if(mb==3),return,end
  if(mb==99),  % the letter "c" on the DEC3100 keyboard
%   copy last cdf profile to new location
    icopy=icopy+1;
    cdfout=cdf;
    cdfout(2:3)=int2str(icopy);
    eval(['! cp ' cdf ' ' cdfout]);
    cdfid2=mexcdf('open',cdfout,'write');
    mexcdf('varput1',cdfid2,'latitude',0,y);
    mexcdf('varput1',cdfid2,'longitude',0,-x);
    mexcdf('close',cdfid2);
    lon(m+icopy-iinit)=x;
    lat(m+icopy-iinit)=y;
    sta(m+icopy-iinit,:)=cdfout;
    hold on
    plot(x,y,'c7x')
    hold off
else
  ind=near(lon+i*lat,x+i*y,1);
%% open the station file and read the temp and salin data
  cdf=sta(ind,:) 
  cdfid=mexcdf('open',cdf,'nowrite');
  [depnam,dlen]=mexcdf('diminq',cdfid,'depth');
  depth=mexcdf('varget',cdfid,'depth',0,dlen);
  corner=[0 0 0 0];
  edge=[1 dlen 1  1];
  u1=mexcdf('varget',cdfid,var1,corner,edge);
  lab1=mexcdf('attget',cdfid,var1,'long_name');
  xlab1=mexcdf('attget',cdfid,var1,'units');
  u2=mexcdf('varget',cdfid,var2,corner,edge);
  lab2=mexcdf('attget',cdfid,var2,'long_name');
  xlab2=mexcdf('attget',cdfid,var2,'units');
  u3=mexcdf('varget',cdfid,var3,corner,edge);
  lab3=mexcdf('attget',cdfid,var3,'long_name');
  xlab3=mexcdf('attget',cdfid,var3,'units');
  rcode=mexcdf('close',cdfid);
  if(first),
   clearg
   subplot(221)
   axis('square')
   first=0;
  end
  plot(u1,-depth)
  ylabel('Depth (m)')
  xlabel(xlab1);
  title(lab1);
  plot(u2,-depth)
  ylabel('Depth (m) ')
  xlabel(xlab2)
  title(lab2)
  plot(u3,-depth)
  ylabel('Depth (m) ')
  xlabel(xlab3)
  title(lab3)
  axis(ax);
  plot(lon,lat,'gx',coast(:,1),coast(:,2),'r-',lon(ind),lat(ind),'wo')
  axis;
  end
end
