% run gen_apex_data first (calls rd_apex), which opens apex edf files, gets
    % data, and saves apex_data.mat
% loads apex_data.mat, opens dat files, and consolidates data from each into array alldat
% applies instrument calibration

% regresses data
% calculates slope for bb calibration value
  
% written by E. Boss and L. Taylor, UMaine, 10/19/06

close all
clear all

load apex_data

for i=1:N % for all files (pos, dat)
    file_i=int2str(i); % change i to string and store as file_i
    n=length(file_i); % check length of file_i
        if n==1 %if 1, 2, 3...9
            name=['pos00',file_i]; % concatenates pos00 and i (pos001) and assigns to array name
            name2=['dat00',file_i]; % concatenates dat00 and i (dat001) and assigns to array name2
            name3=['date00',file_i];
        elseif n==2 %if 10, 11, 12...99 
             name=['pos0',file_i]; % concatenates pos0 and i (pos010) and assigns to array name
             name2=['dat0',file_i]; % concatenates dat0 and i (dat010) and assigns to array name2
             name3=['date0',file_i];
        else %if 100, 101, 102...999
             name=['pos',file_i]; % concatenates pos and i (pos100) and assigns to array name
             name2=['dat',file_i]; % concatenates dat and i (dat100) and assigns to array name2
             name3=['date',file_i];
        end           
    A=eval(name); % eval(pos001) - takes strings in pos001 (lat, lon, minutes) and changes to numbers, then assigns to array A
    latlong=A(1,1:2); 
    long=[latlong(:,1)-60]; % longitude
    lat=latlong(:,2); % latitude
    B=eval(name2); % data
    C=eval(name3); % date
    [mm nn]=size(B);

    % open each file (pos, dat and date) 000, 001, 000+ in turn
    if i==1 % if first time through (pos001, dat001 and date001)
        dat=[B,zeros(mm,1)]; % create array dat from data in dat file with last column (16) of zeros
        alldat=[dat,zeros(mm,1:2)]; % create array alldat with dat plus columns 17 and 18 of zeros
        alldat(:,17)=long(1,1); % replace column 17 with longitude
        alldat(:,18)=lat(1,1); % replace column 18 with latitude
         
    else %if not the first time around
        AA=[B,(i-1)*5*ones(mm,1)]; % create array AA from file and replace/increment last column (16) to get profile stages 
        dat=[dat;AA]; % append AA to dat
        CC=[AA,zeros(mm,1:2)]; % create array CC from AA, and add two columns of zeros
        CC(:,17)=long(1,1); % replace column 17 with longitude
        CC(:,18)=lat(1,1); % replace column 18 with latitude
        alldat=[alldat;CC]; % append CC to alldat array
    end
end

% remove negative values (no data collected)
I=find(alldat<-100);
alldat(I)=NaN;

%variables are:
P=alldat(:,1); % depth (dB)
T=alldat(:,2); % temperature (degrees C)
S=alldat(:,3); % salinity (psu)
O=alldat(:,4);
BlueRef=alldat(:,5);
BlueSig=alldat(:,6);
LSSRef=alldat(:,7);
LSSSig=alldat(:,8);
theta=alldat(:,9); % temperature (degrees C)
sigma=alldat(:,10); % density-1000 (kg/m3)
O2sat=alldat(:,11); % oxygen (% saturation)
O2=alldat(:,12); % oxygen (um/kg)
Chl=alldat(:,13); % chlorophyll (mg/m3, not calibrated)
B_LSS=alldat(:,14); % bb (m-1, not calibrated)
date=alldat(:,15); % date (datenum format)
profile=alldat(:,16); % profile (in days from start)
long=alldat(:,17); % longitude (degrees)
lat=alldat(:,18); % latitude (degrees)

% calibrations:
B_LSS_b=(LSSSig-51.8)*0.00004637; % instrument calibration
cp=B_LSS*63;
cp_b=B_LSS_b*63;
Chl_b=(BlueSig-25)*0.01079772; % instrument calibration; 25 is mean of counts at >950m

% get profile number
profile=profile/5+1; 

% get time beginning June 17, 2004
tt=datenum(2004,6,17)+alldat(:,16);
t=168+alldat(:,16);

% make array of all bb data
calBb=[date,P,T,S,B_LSS_b,tt];

% get bb data from 1 to 9m
I=find(calBb(:,2)>1 & calBb(:,2)<9);
floatDate=calBb(I,1); % 216x1
floatT=calBb(I,3); % 216x1
floatS=calBb(I,4); % 216x1
floatBb=calBb(I,5); % 216x1
tt=tt(I); % 216x1

% clear unnecessary variables from workspace
clear dat dat0* dat1* dat2* date0* date1* date2* pos* alldat;
clear name* filename* file*;
clear A* i* B C CC N mm n* t;
clear latlong lat long;
clear B_LSS B_LSS_b BlueRef BlueSig Chl Chl_b Chl_b_corr LSSRef LSSSig O O2 O2sat P POC_b C_phy;
clear cp cp_b profile sigma theta;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIS data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\0005\MODIS\MODIS_spectra_7_5km_dilate.mat';
% MODIS_spectra_7_5km_dilate.mat contains the following variables within 7.5km of the float:
% wl
% mean_chl
% median_chl
% std_chl
% date_modis
% spectra_median
% spectra_mean
% spectra_std

Lwn=spectra_median;
F0=[172.718,187.688,194.973,185.820,186.594];
mm=length(Lwn(:,1));
Rrs=1.38*Lwn(:,1:5)./(ones(mm,1)*F0);
II=find(isfinite(Rrs(:,1).*Rrs(:,2).*Rrs(:,3).*Rrs(:,4).*Rrs(:,5)));
Rrs=Rrs(II,:);
time_MODIS=date_modis(II); % 265x1
wavelength=[412 443 488 531 551];

%user inputs:
%Temperature and salinity of water under consideration [degree C] (to compute water IOPs)
Temp=interp1(tt,floatT,time_MODIS,'linear'); % 265x1
Salt=interp1(tt,floatS,time_MODIS,'linear'); % 265x1

for i=1:length(II)
    %compute IOPs of seawater
    [a_w,bb_w]=IOP_of_sw(wavelength(1:5),Temp(i),Salt(i));
    %a_w=[0.02027,0.0111,0.01677,0.04496,0.05492];%,0.4920];
    %bb_w=[0.003242,0.002454,0.001622,0.001132,0.000989];%,0.0004262];
    rrs=Rrs./(0.52+1.7*Rrs);
    A=GSM01_invert(rrs(i,1:5),a_w,bb_w);
    chl(i)=A(1);
    a_dg(i)=A(2);
    b_bp(i)=A(3);
end

figure(1); % wavelength 443
plot(time_MODIS,b_bp,'g*')
set(gca,'XTick',[732099 732464 732829 733194 733560]);
datetick(gca,'x','mmmyy','keepticks');
ylabel('b bp (443)');

%%%%%%%%%%%%%%%%% CORRELATIONS
%%%%%% median MODIS chl vs float chl
float_sat_median_bb=interp1(time_MODIS, b_bp, floatDate, 'linear'); % size: 216

I=finite(float_sat_median_bb); % size: 203
A1_sat=float_sat_median_bb(I);
B1_float=floatBb(I);

I2=find(A1_sat<0.06); % remove one outlier size: 202
A2=corrcoef(A1_sat(I2), B1_float(I2)); % .9017

%%%% check plot
% loglog(float_sat_median_bb, floatBb,'.');
scatter(A1_sat(I2), B1_float(I2),'.');

%%%%%%%%%%%%%%%%% Regression analysis
sat_bb=A1_sat(I2);
float_bb=B1_float(I2);

[poly_bb1] = polyfit(sat_bb(~isnan(sat_bb)), ...
                    float_bb(~isnan(sat_bb),1), 1);

[poly_bb2] = polyfit(float_bb(~isnan(sat_bb),1), ...
                     sat_bb(~isnan(sat_bb)),1);

slope_bb=sqrt(poly_bb1(1)*(1/poly_bb2(1))); % 2.8127 (1.8903 dilated) USE THIS VALUE FOR CALIBRATION

[poly_bb3] = polyfit(log(sat_bb(~isnan(sat_bb))), ...
                    log(float_bb(~isnan(sat_bb),1)), 1);

[poly_bb4] = polyfit(log(float_bb(~isnan(sat_bb),1)), ...
                     log(sat_bb(~isnan(sat_bb))),1);

slope_bb_2=exp((poly_bb3(2)-poly_bb4(2))/2); % (1.5259 dilated) DO NOT USE THIS VALUE 

plot(floatDate(~isnan(sat_bb)),sat_bb(~isnan(sat_bb)),'*',floatDate(~isnan(sat_bb)),float_bb(~isnan(sat_bb),1)/2.8127,'o');
