%% EXPORTS
clear all
close all
filename='partE3_6.mat';
load(['/Users/emmanuel.boss/Desktop/EXPORTS01/IOPPackage/software/Particulate_output_E1_1_7_E2_1_8_E3_1_6/',filename])
%output file
textfile1='EXPORTS-EXPORTSNP_SlowDROP-ACS_process_20180907b_ap_cp_R1.sb';
textfile2='EXPORTS-EXPORTSNP_SlowDROP-ACS_process_20180907b_Products_R2.sb';
%This section insures data to be output is reasonable
% ref = find(wl>=440, 1,'first');
% I=find(ap(:,ref)<0.01);
% ap(I,:)=NaN;
ref = find(wl>=734, 1,'first');
I=find(abs(ap(:,ref))>0.003);
ap(I,:)=NaN;

ref = find(wl>=650, 1,'first');

% Derive Chl (Line heigh at 676 compared to 650 and 715)
% 
ap_lh = interp1(wl, ap',[650 676 715],'linear')';
line_height = (ap_lh(:,2)-(39/65*ap_lh(:,1)+26/65*ap_lh(:,3)));
chl_lh = 138.14 * line_height.^1.11; % EXPORTS relation

%   The particulate organic carbon (POC) is computed using the particulate attenuation at 660 nm Using the global relationship from Gardner et al. (2006)
cp660 = interp1(wl,cp',660,'linear')';
cp_poc = cp660.*380;

figure(1)
subplot(2,1,1)
plot(wl,cp);
ylabel('c_p(\lambda)[m^-1]')
xlabel('wavelength [nm]')
subplot(2,1,2)
plot(wl,ap);
ylabel('a_p(\lambda)[m^-1]')
xlabel('wavelength [nm]')

export_fig([filename(1:8),'a.png'])

opts = optimset('fminsearch');
%opts = optimset(opts,'NonlEqnAlgorithm', 'gn');
opts = optimset(opts,'MaxIter',20000000);
opts = optimset(opts,'MaxFunEvals',20000);
opts = optimset(opts,'TolX',1e-8);
opts = optimset(opts,'TolFun',1e-8);
II = find(wl<=700);
ref2 = find(wl>=650, 1,'first');
x1=ones(length(P),2)*NaN;
for i=1:length(P)
    if isfinite(cp(i,ref2))
        [x1(i,:), fiterr(i)] = fminsearch(@least_squares_cp, [cp(i,ref2) 1], opts, cp(i,II), wl(II));
    end
end
figure(2)
cp_gamma=x1(:,2);
subplot(2,2,1)
plot(cp_poc,-P,'*');
ylabel('depth')
xlabel('cp_{poc}')
subplot(2,2,2)
plot(chl_lh,-P,'*');
ylabel('depth')
xlabel('chl_{lh}')
subplot(2,2,3)
plot(cp_gamma,-P,'*');
ylabel('depth')
xlabel('\gamma')
export_fig([filename(1:8),'b.png'])

II=find(wl>440 & wl < 705);
amps=ones(length(P),13)*NaN;
for i=1:length(P)
    if isfinite(ap(i,20))
        I=find(isfinite(ap(i,:)));
        if length(I)>40
             [amps(i,:),compspec,sumspec] = spectral_decomp(wl(I),ap(i,I),wl(I),unc_ap(I,i),1);
             err(i)=sum(abs(ap(i,II)-sumspec(II)'))/length(II);
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%inputs for sb file:

docfile1='EXPORTS_ap_cp_processing_SlowDROP.docx';

roundwl=roundn(wl,-1);
for iwl=1:length(wl)
    wlarray{iwl}=['ap',num2str(roundwl(iwl)),','];
    wlarray2{iwl}=['ap',num2str(roundwl(iwl)),'_sd,'];
    wlarray3{iwl}=['cp',num2str(roundwl(iwl)),','];
end
for iwl=1:length(wl)-1
    wlarray4{iwl}=['cp',num2str(roundwl(iwl)),'_sd,'];
end
wlarray4{length(wl)}=['cp',num2str(roundwl(length(wl))),'_sd'];

% create arrays of units for the file header
apunits=repmat('1/m,',1,4*length(wl)-1);
apunits=[apunits,'1/m'];

%wlprint = strrep(num2str(wlem'),'       ',' ')

fileID(1)=fopen(textfile1,'w');
fileID(2)=fopen(textfile2,'w');

%header:
for i=1:2
fprintf(fileID(i),'%-30s\n','/begin_header');
fprintf(fileID(i),'%-30s\n','/investigators=Emmanuel_Boss,Nils_Haentjens');
fprintf(fileID(i),'%-30s\n','/affiliations=UMaine');
fprintf(fileID(i),'%-30s\n','/contact=emmanuel.boss@maine.edu');
fprintf(fileID(i),'%-30s\n','/experiment=EXPORTS');
fprintf(fileID(i),'%-30s\n','/cruise=EXPORTSNP');
fprintf(fileID(i),'%-30s\n','/station=NA');
if i==1
   fprintf(fileID(i),'%-30s\n',['/data_file_name=',textfile1]);
else
   fprintf(fileID(i),'%-30s\n',['/data_file_name=',textfile2]);
end
fprintf(fileID(i),'%-30s\n','/documents=EXPORTS_ap_cp_processing_SlowDROP.docx');
fprintf(fileID(i),'%-30s\n','/calibration_files=acs298.dev');
fprintf(fileID(i),'%-30s\n','/data_type=cast');
% fprintf(fileID(i),'%-30s\n','/data_status=final');
fprintf(fileID(i),'%-30s\n','/start_date=20180907');
fprintf(fileID(i),'%-30s\n','/end_date=20180907');
fprintf(fileID(i),'%-30s\n','/start_time=19:53:00[GMT]');
fprintf(fileID(i),'%-30s\n','/end_time=20:17:00[GMT]');
fprintf(fileID(i),'%-30s\n','/north_latitude=50.413[DEG]');																														
fprintf(fileID(i),'%-30s\n','/south_latitude=50.4125[DEG]');
fprintf(fileID(i),'%-30s\n','/east_longitude=-144.6729[DEG]');
fprintf(fileID(i),'%-30s\n','/west_longitude=-144.6763[DEG]');
fprintf(fileID(i),'%-30s\n','/water_depth=-999');
fprintf(fileID(i),'%-30s\n','/delimiter=comma');
fprintf(fileID(i),'%-30s\n','/missing=-999');
fprintf(fileID(i),'%-30s\n','/instrument_model=ac-s');
fprintf(fileID(i),'%-30s\n','/instrument_manufacturer=Sea-Bird_Scientific');
end
%spectral output
fprintf(fileID(1),'%-30s\n','!COMMENTS');
fprintf(fileID(1),'%-30s\n','!ap and cp values are based on binning 1m measurements; calculations explained in associated documents');
fprintf(fileID(1),'%-300s\n',['/fields=pressure,Wt,sal,',wlarray{:},wlarray2{:},wlarray3{:},wlarray4{:}]);
fprintf(fileID(1),'%-300s\n',['/units=dbar,degreesC,PSU,',apunits]);
fprintf(fileID(1),'%-30s\n','/end_header');

%product output
fprintf(fileID(2),'%-30s\n','!COMMENTS');
fprintf(fileID(2),'%-30s\n','!see associated documents for products');
fprintf(fileID(2),'%-300s\n',['/fields=pressure,Wt,sal,Chl_lineheight,POC_cp,cp_gamma']);
fprintf(fileID(2),'%-300s\n',['/units=dbar,degreesC,PSU,ug/L,ug/L,unitless']);
fprintf(fileID(2),'%-30s\n','/end_header');

%define the format for the data part formdate=repmat('%9s ',1,3);
formrad=repmat('%6.4f,',1,(length(wl)*4)-1);
format=['%6.4f,','%6.4f,','%6.4f,',formrad,'%6.4f\n'];

%delete lines with no data
A=[ap,cp,unc_ap',unc_cp'];
B=A*NaN;
index=[];
k=1;
for i=1:length(P)
    if isfinite(A(i,40)) | isfinite(A(i,120))
        index(k)=i;
        I=find(isnan(A(i,:)));
        A(i,I)=-999;
        k=k+1;
    end
end
I=isfinite(cp_poc)==0;
cp_poc(I)=-999;
I=isfinite(chl_lh)==0;
chl_lh(I)=-999;
I=isfinite(cp_gamma)==0;
cp_gamma(I)=-999;

data1=[P(index)',T(index),S(index),A(index,:)];
data2=[P(index)',T(index),S(index),chl_lh(index),cp_poc(index),cp_gamma(index)];

format2=['%6.4f,','%6.4f,','%6.4f,','%6.4f,','%6.4f,','%6.4f\n'];

for ii = 1:length(index)
    fprintf(fileID(1),format,data1(ii,:));
end

fclose(fileID(1))

for ii = 1:length(index)
    fprintf(fileID(2),format2,data2(ii,:));
end

fclose(fileID(2))

