% convert CTD data from structure to matrix format
warning off all
clear all;

stations=[2 4:39 41:62];

for zz=1:length(stations)
    ss=stations(zz)

%load in each .cnv processed file from each cast
%put needed data into a structure
if ss<10
    
    [a1,lat_deg,lat_min,ns]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_00'... 
        int2str(ss) '.cnv'],...
        '%18c%2n%5f%c', 1, 'headerlines',9);
    
    [a1,lon_deg,lon_min,ew]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_00'... 
        int2str(ss) '.cnv'],...
        '%19c%3n%5f%c', 1, 'headerlines',10);
    
    [a1,month,month_day,year,hour,a2,minute,a3,sec]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_00'... 
        int2str(ss) '.cnv'],...
        '%20c%s%n%n%n%c%n%c%n', 1, 'headerlines',11);
    
    [time,pres,temp1,temp2,cond1,cond2,oxy1,oxy2,altm,dep,dens1,dens2,theta1,theta2,sound1,sound2,...
        oxy_ml1,oxy_mg1,oxy_mol1,oxy_ml2,oxy_mg2,oxy_mol2,sal1,sal2,flag]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_00'... 
        int2str(ss) '.cnv'],...
        '%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n', 'headerlines',366);
else
    
    [a1,lat_deg,lat_min,ns]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_0'... 
        int2str(ss) '.cnv'],...
        '%18c%2n%5f%c', 1, 'headerlines',9);
    
    [a1,lon_deg,lon_min,ew]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_0'... 
        int2str(ss) '.cnv'],...
        '%19c%3n%5f%c', 1, 'headerlines',10);
    
    [a1,month,month_day,year,hour,a2,minute,a3,sec]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_0'... 
        int2str(ss) '.cnv'],...
        '%20c%s%n%n%n%c%n%c%n', 1, 'headerlines',11);
    
    [time,pres,temp1,temp2,cond1,cond2,oxy1,oxy2,altm,dep,dens1,dens2,theta1,theta2,sound1,sound2,...
        oxy_ml1,oxy_mg1,oxy_mol1,oxy_ml2,oxy_mg2,oxy_mol2,sal1,sal2,flag]=...
        textread(['/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/ctd/pne17/hydro/proc_data/dpne17_0'... 
        int2str(ss) '.cnv'],...
        '%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n%n', 'headerlines',366);
end
lat=lat_deg+(lat_min/60);
if ns=='S'
    lat=-lat;
end
lon=lon_deg+(lon_min/60);
if ew=='W'
    lon=-lon;
end

if sum(month{1}=='Jan')==3
    month=1;
elseif sum(month{1}=='Feb')==3
    month=2;
elseif sum(month{1}=='Mar')==3
    month=3;
elseif sum(month{1}=='Apr')==3
    month=4;
elseif sum(month{1}=='May')==3
    month=5;
elseif sum(month{1}=='Jun')==3
    month=6;
elseif sum(month{1}=='Jul')==3
    month=7;
elseif sum(month{1}=='Aug')==3
    month=8;
elseif sum(month{1}=='Sep')==3
    month=9;
elseif sum(month{1}=='Oct')==3
    month=10;
elseif sum(month{1}=='Nov')==3
    month=11;
elseif sum(month{1}=='Dec')==3
    month=12;
end

%disp(max(pres));
cast_num=ss*ones(length(pres),1);

cnv_dlp=struct('pressure',pres,'temppri',temp1,'oxmllpri',oxy_ml1,'salpri',sal1,'station',cast_num,...
    'latitude',lat,'longitude',lon,'year',year,'month',month,'month_day',month_day,'hour',hour,...
    'minute',minute','sec',sec);

% test for bad salinity %
for si=1:length(pres)
    if (cond1(si)<=0),
        disp('BAD CONDUCTIVITY');
        disp(zz);
        disp(si);
%         cnv_dlp.salpri(1:si)=nan;
%         cnv_dlp.temppri(1:si)=nan;
%         cnv_dlp.oxmllpri(1:si)=nan;        
        cnv_dlp.salpri(si)=nan;
        cnv_dlp.temppri(si)=nan;
        cnv_dlp.oxmllpri(si)=nan;
    end
end

% %test for static stability - Sometimes it fails this test, but I think
% this should be fixed after the fact.
% for si=1:length(pres)-1
%     if (dens1(si)>dens1(si+1) & pres(si)<500)
%         disp('UNSTABLE');
%         disp(zz);
%         disp(si);
%         cnv_dlp.salpri(si)=nan;
%         cnv_dlp.temppri(si)=nan;
%         cnv_dlp.oxmllpri(si)=nan;
%     end
% end


if (1)  % apply calibration coefficents derived by
        % calibrate_sal_oxygen.m.  The model for the
        % bottle/sensor differences is hard-coded here. 
        
    %disp('Applying calibration adjustment to primary sal, O2 sensor profiles') 
    load calibration_coeff xs xo;
    cnv_cal=cnv_dlp;
    %for j=1:length(cnv_cal);
    for j=1:1 
        cnv_cal(j).salpri=cnv_cal(j).salpri - ( ...
            xs(1)+xs(2)*cnv_cal(j).temppri + xs(3)*(cnv_cal(j).temppri).^2 + ...
            xs(4)*cnv_cal(j).pressure + xs(5)*(cnv_cal(j).pressure).^2 + ...
            xs(6)*cnv_cal(j).station + xs(7)*(cnv_cal(j).station).^2 + ...
            xs(8)*cnv_cal(j).salpri + xs(9)*(cnv_cal(j).salpri).^2  );
        
        cnv_cal(j).oxmllpri=cnv_cal(j).oxmllpri - (...
                xo(1)+xo(2)*cnv_cal(j).temppri + xo(3)*(cnv_cal(j).temppri).^2 + ...
                xo(4)*cnv_cal(j).pressure + xo(5)*(cnv_cal(j).pressure).^2 + ...
                xo(6)*cnv_cal(j).station + xo(7)*(cnv_cal(j).station).^2 + ...
                xo(8)*cnv_cal(j).oxmllpri + xo(9)*(cnv_cal(j).oxmllpri).^2  );
    end;
    ctd(zz)=cnv_cal;
    calflag=1;
else    
    eval(['ctd=',sname,';']);
    calflag=0;
end;

% NOTE: The calibration for equatorial casts isn't good below 1500m, the actual
% values are not far from the bottle values and the model over corrected
% the values!
% Slope comes from diff(samples(q(1:2),11)-samples(q(1:2),5))/(samples(q(1),9)-samples(q(2),9)) - ((0.1205/2)/(samples(q(1),9)-samples(q(2),9)))
if ss==16,
%    ctd(zz).oxmllpri(1993:3492)=oxy_ml1(1993:3492)-4.8212e-05*(pres(1993:3492)-pres(1993))+0.1205/2;  % Try to linearly correct the pressures.
 ctd(zz).oxmllpri(1900:3492)=oxy_ml1(1900:3492)-7.9583e-06*(pres(1900:3492)-pres(1900))+0.1205/2;  % Try to linearly correct the pressures.
end;

%clear sname;

% generate plot of uncalibrated sensor, subsampled sensor values
% at bottle depths, bottle values, and calibrated salinity/oxygen 
%load calibrations samples

% load('/Users/renellys.perez/Breve/Observations/PNE2017_RB1701/calibration/results/calibration_data.mat');
% 
% stn=zz;
% close; 
% subplot(121)
% plot(sal1,-pres,'r','LineWidth',1);
% hold on;
% plot(ctd(stn).salpri,-ctd(stn).pressure,'b','LineWidth',2);
% q=find(samples(:,1)==ctd(stn).station(1));
% % plot sensor values as red dots (red dots should fall exactly on red curve)
% plot(samples(q,3),-samples(q,9),'.r');
% % plot autosal values as blue circles (blue curve should generally follow)
% plot(samples(q,10),-samples(q,9),'ob');
% hold off;
% title(['Cast ',num2str(ctd(stn).station(1)),' salinity: red uncal, blue cal'])
% 
% subplot(122)
% plot(oxy_ml1,-pres,'r','LineWidth',1);
% hold on;
% plot(ctd(stn).oxmllpri,-ctd(stn).pressure,'b','LineWidth',2);
% q=find(samples(:,1)==ctd(stn).station(1));
% % plot sensor values as red dots (red dots should fall exactly on red curve)
% plot(samples(q,5),-samples(q,9),'.r');
% % plot titration values as blue circles (blue curve should generally follow)
% plot(samples(q,11),-samples(q,9),'ob');
% hold off;
% title(['Cast ',num2str(ctd(stn).station(1)),' oxygen: red uncal, blue cal'])
% %end;
% 
% printstring=['print -depsc /Users/renellys.perez/Breve/Observations/PNE2017_RB1701/calibration/results/check_calibrated_curves_', int2str(ss) '.eps']
% eval(printstring);

end     %loop through each cast



%stnnbr=(1:length(ctd))';
stnnbr=stations;
%stnnbr=stnnbr-1; % cast 0 = first entry

p=(0:3753)';        %max pressure from all casts

temp=NaN*ones(length(p),length(stnnbr)-1);
sal=temp; oxy=temp;
lat=NaN*ones(length(stnnbr)-1,1);
lon=lat;  t=lon;

for zz=1:length(stations);
    stn=stations(zz);
    q=find(stnnbr==stn);
    lat(q)=ctd(q).latitude;
    lon(q)=ctd(q).longitude;
    t(q)=julian([ctd(q).year ctd(q).month ctd(q).month_day ctd(q).hour ctd(q).minute ctd(q).sec]);
    ttemp=ctd(q).temppri;
    tp=ctd(q).pressure;
    temp(:,q)=interp1(tp,ttemp,p','linear')';
    tsal=ctd(q).salpri;
    sal(:,q)=interp1(tp',tsal,p','linear')';
    toxy=ctd(q).oxmllpri;
    oxy(:,q)=interp1(tp',toxy,p','linear')';
end; 
clear stn q ttemp tsal toxy;

% Define coordinates

xdesc='Latitude';
ydesc='Pressure (db)';
xrange=[-5 20.5]
yrange=[-1500 0];
xtraline=[0 0];
ytraline=[yrange]; xtratype='k--'; xtrawidth=2;
axis_scale=0.2;
xinches=5;
yinches=8*axis_scale;
ticksx=[-5:2.5:20];
ticksy=[-1600:200:0];
ticksx_label=['5S||EQ||5N||10N||15N||20N|'];
ticksy_label=['1600||1200||800||400||0'];

caxis1=[34 37];
cint1=[30 34:0.25:37 3]; cint1bar=[34:0.25:37];
cint1b=[34:0.25:37];
cint1c=[34:0.25:37];
bold1=[0 0];
caxislab=['S (psu)'];
cxtick1=[34:37];
cxlabel1=['34|35|36|37'];


close all;
figure('PaperUnits','inches','Units','inches','Position',[2 2 11 8]);
set(gcf,'PaperPosition',[0.25 0.25 10.5 8]);
h1=axes('Units','Inches','Position',[1.5 3.25 xinches yinches]);
contourf(lat(1:58),-p,sal(:,1:58),cint1)
caxis(caxis1);
hold on;
axis([-5 20.5 -1500 0]);
title('PSAL','Fontsize',16,'Fontname','Palatino');
plot(xtraline,ytraline,xtratype,'LineWidth',xtrawidth);
ylabel(ydesc,'Fontsize',14,'Fontname','Palatino');
xlabel('','Fontsize',14,'Fontname','Palatino');
axis([xrange yrange])
set(h1,'XTick',ticksx,'XTickLabel',ticksx_label,'Fontsize',14,'Fontname','Palatino');
set(h1,'YTick',ticksy,'YTickLabel',ticksy_label,'Fontsize',14,'Fontname','Palatino','TickDir','out');

% COLORBAR:

hcol=axes('Units','Inches','Position',[6.75 3.25 0.15 yinches]);
pcolor([-1 1],cint1bar,([1 -1]'*cint1bar)');
caxis(caxis1);
set(hcol,'XTick',[],'XTickLabel',[],'Fontsize',10,'Fontname','Palatino');
set(hcol,'YTick',[cxtick1],'YTickLabel',cxlabel1,'Fontsize',12,'Fontname','Palatino','YAxisLocation','right');
ylabel(caxislab,'Fontsize',12,'Fontname','Palatino');



caxis1=[0 6];
cint1=[0:0.3:6]; cint1bar=[0:0.3:6];
cint1b=[0:0.3:6];
cint1c=[0:0.3:6];
bold1=[0 0];
caxislab=['O_2 (mL/L)'];
cxtick1=[0:2:6];
cxlabel1=['0|2|4|6'];


h1=axes('Units','Inches','Position',[1.5 1 xinches yinches]);
contourf(lat(1:58),-p,oxy(:,1:58),[0:0.3:6])
caxis(caxis1);
hold on;
axis([-5 20.5 -1500 0]);
title('Dissolved Oxygen','Fontsize',16,'Fontname','Palatino');
plot(xtraline,ytraline,xtratype,'LineWidth',xtrawidth);
ylabel(ydesc,'Fontsize',14,'Fontname','Palatino');
xlabel('','Fontsize',14,'Fontname','Palatino');
axis([xrange yrange])
set(h1,'XTick',ticksx,'XTickLabel',ticksx_label,'Fontsize',14,'Fontname','Palatino');
set(h1,'YTick',ticksy,'YTickLabel',ticksy_label,'Fontsize',14,'Fontname','Palatino','TickDir','out');

% COLORBAR:

hcol=axes('Units','Inches','Position',[6.75 1 0.15 yinches]);
pcolor([-1 1],cint1bar,([1 -1]'*cint1bar)');
caxis(caxis1);
set(hcol,'XTick',[],'XTickLabel',[],'Fontsize',10,'Fontname','Palatino');
set(hcol,'YTick',[cxtick1],'YTickLabel',cxlabel1,'Fontsize',12,'Fontname','Palatino','YAxisLocation','right');
ylabel(caxislab,'Fontsize',12,'Fontname','Palatino');
% 
% % Color
% 
% map=lbmap(20,'BrownBlue');
% colormap(flipud(map));

print -depsc pne2017_ctd.eps
    
save ctddata_pne2017 lon lat t p temp sal oxy calflag stnnbr;

