% run gen_msg_data first (calls rd_msg) to open msg files and get data
% this file loads apex_msg file and consolidates data from individual arrays into one array (dat)
% adds stat, lon, lat, date, and date2 data to array dat
% applies instrument calibration
% gets first float reading of every day
% loads modis data and gets median of chl readings within 7.5km of float
% regresses data (linear)
% calculates slope for calibration value

clear all
close('all')

load apex_msg % from gen_msg_data
for i=1:N %for files 1 through end - N from gen_msg_data
    file_i=int2str(i); %change i to string and store as file_i
    n=length(file_i); %gets length of string file_i
        if n==1 %if length is 1 (i is 1,2,3....9)
            name2=['dat00',file_i]; %concatenate dat00 and file_i (dat001, dat002, dat003...dat009) and store in array name2
        elseif n==2 %if length is 2 (i is 10, 11, 12...99)
             name2=['dat0',file_i]; %concatenate dat0 and file_i (dat010, dat011, dat012...dat099) and store in array name2
        else %if length is 2 (i is 100, 101, 102...999)
             name2=['dat',file_i]; %concatenate dat and file_i (dat100, dat101, dat102...dat999) and store in array name2
        end
    B=eval(name2); %takes strings in dat001 and changes to numbers, then assigns to array B
    [mm nn]=size(B); %returns number of rows and columns in B (dat001)
    if i==1 %if this is the first time through
        dat=[B,zeros(mm,5)]; %make array dat with B and five extra columns of zeros (for as many rows as there are in B)
    else %if this is not the first time through
        AA=[B,zeros(mm,5)]; %make array AA with B and five extra columns of zeros
        dat=[dat;AA]; %append array AA to dat and make array AA - should have 6 data fields plus empty columns
    end
end
%stat combine stat files into one
for i=1:N
    file_i=int2str(i);
    n=length(file_i);
        if n==1
             name2=['stat00',file_i];
        elseif n==2
             name2=['stat0',file_i];
        else
             name2=['stat',file_i];
        end
    B=eval(name2);
    [mm nn]=size(B);
    if i==1
        stat=[B];
    else
       AA=[B];
       stat=[stat;AA];
    end
end
%lat combine lat files into one
for i=1:N
    file_i=int2str(i);
    n=length(file_i);
        if n==1
             name2=['lat00',file_i];
        elseif n==2
             name2=['lat0',file_i];
        else
             name2=['lat',file_i];
        end
    B=eval(name2);
    [mm nn]=size(B);
    if i==1
        lat=[B];
    else
       AA=[B];
       lat=[lat;AA];
    end
end
%lon combine lon files into one
for i=1:N
    file_i=int2str(i);
    n=length(file_i);
        if n==1
             name2=['lon00',file_i];
        elseif n==2
             name2=['lon0',file_i];
        else
             name2=['lon',file_i];
        end
    B=eval(name2);
    [mm nn]=size(B);
    if i==1
        lon=[B];
    else
       AA=[B];
       lon=[lon;AA];
    end
end
%date combine date files into one
for i=1:N
    file_i=int2str(i);
    n=length(file_i);
        if n==1
             name2=['date00',file_i];
        elseif n==2
             name2=['date0',file_i];
        else
             name2=['date',file_i];
        end
    B=eval(name2);
    [mm nn]=size(B);
    if i==1
        date=[B];
    else
       AA=[B];
       date=[date;AA];
    end
end
%date2 combine date2 files into one
for i=1:N
    file_i=int2str(i);
    n=length(file_i);
        if n==1
             name2=['date200',file_i];
        elseif n==2
             name2=['date20',file_i];
        else
             name2=['date2',file_i];
        end
    B=eval(name2);
    [mm nn]=size(B);
    if i==1
        date2=[B];
    else
       AA=[B];
       date2=[date2;AA];
    end
end
%
%add stat, lat, lon, date and date2 to array dat
dat(:,6)=str2num(stat); %replaces column 6 zeros with stat
dat(:,7)=str2num(lat); %replaces column 7 zeros with lat
dat(:,8)=str2num(lon);  %replaces column 8 zeros with lon
dat(:,9)=date;  %replaces column 9 with zeros date
dat(:,10)=date2;  %replaces column 10 with zeros date2

P=dat(:,1); %get pressure from array dat(column 1)and store in P
FRef=dat(:,2); %F_Ref
BlueSig=dat(:,3); %F_Sig (blue sig = chlorophyll)
LSSRef=dat(:,4); %LSS_Ref
LSSSig=dat(:,5); %LSS_Sig (backscattering)
stat=dat(:,6); % stat
lat=dat(:,7); % lat
lon=dat(:,8); % lon
date=dat(:,9); % date
date2=dat(:,10); % date2

% calibrations:
Chl_b=(BlueSig-25)*0.01079772; % instrument calibration; 25 is mean of counts at >950m

% make array of all surface data
calChl=[date,Chl_b]; 

%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTS
% plot chl to check
% figure(1);
% plot(date,Chl_b,'.')
% datetick('x','mmmyy')
% title({'Chlorophyll','Float 0005, North Atlantic'})
% ylabel('Concentration [mg/m^3]')
% xlabel('time [month]')
%%%%%%%%%%%%%%%%%%%%%%%%%%

% clear unnecessary variables from workspace
clear dat dat0* dat1* dat2* date0* date* pos*;
clear AA B B_LSS_b B_LSS_b_corr BlueSig C Chl_b Chl_b_corr D E F FRef G H LSSRef LSSSig N P;
clear file_i file_list_out filename filename1 i ii;
clear lat lat0* lat1* lon lon0* lon1* lat2* lon2*;
clear mm n n_file_out name2 nn;
clear stat stat0* stat1* stat2*;

% get first reading of the day
I=find(diff(calChl(:,1))>4);
time_surf=[calChl(1,1);calChl(I+1,1)]; % size: 221
chl_surf=[calChl(1,2);calChl(I+1,2)]; % size: 221

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MODIS data
% get modis data
load 'C:\Documents and Settings\Lisa Taylor\My Documents\Data\0005\MODIS\modis_chl_dilate.mat'
% modis_chl_dilate contains:
  % modis_chl_ts (24646x2601) - chl data only from 06/13/04 - 01/01/07 THIS IS A 51X51 GRID
  % modisChl [date, chl] (24646x2602)
  % distgrid (24646x2601) grid of distance around float
modis_date=modisChl(:,1);

%%% get mean and median modis values within 7.5km of float
[uu vv]=size(modis_chl_ts); % 2116x2601
I7km_modis=find(distgrid(1,:)<7.5); 
for i=1:uu
   dist_7km_A_modis(i)=nanmean(modis_chl_ts(i,I7km_modis),2);
   dist_7km_B_modis(i)=nanmedian(modis_chl_ts(i,I7km_modis),2);
end

dist_7km_mean_modis=dist_7km_A_modis';
dist_7km_median_modis=dist_7km_B_modis';

% remove zero's
I7k_modis=find(dist_7km_mean_modis>0); % size: 301
I7kk_modis=find(dist_7km_median_modis>0); % size: 301

%%%%%%%%%%%%%%%%% CORRELATIONS -208pts (missing 13)
%%%%%% median MODIS chl vs float chl
float_sat_median_chl=interp1(modis_date(I7kk_modis), dist_7km_median_modis(I7kk_modis), time_surf, 'linear');

I=finite(float_sat_median_chl);
A1=corrcoef(float_sat_median_chl(I), chl_surf(I)); % .8852
%log correlation is 0.712

%%%% check plot
% loglog(float_sat_median_chl, chl_surf,'.');
% scatter(float_sat_median_chl, chl_surf,'.');

%%%%%%%%%%%%%%%%% Regression analysis
[poly_chl1] = polyfit(float_sat_median_chl(~isnan(float_sat_median_chl)), ...
                    chl_surf(~isnan(float_sat_median_chl),1), 1);

[poly_chl2] = polyfit(chl_surf(~isnan(float_sat_median_chl),1), ...
                     float_sat_median_chl(~isnan(float_sat_median_chl)),1);

slope_chl=sqrt(poly_chl1(1)*(1/poly_chl2(1))); % 5.0359 (4.9281 dilated) USE THIS VALUE FOR CALIBRATION

[poly_chl3] = polyfit(log(float_sat_median_chl(~isnan(float_sat_median_chl))), ...
                    log(chl_surf(~isnan(float_sat_median_chl),1)), 1);

[poly_chl4] = polyfit(log(chl_surf(~isnan(float_sat_median_chl),1)), ...
                     log(float_sat_median_chl(~isnan(float_sat_median_chl))),1);

slope_chl_2=exp((poly_chl3(2)-poly_chl4(2))/2); % 3.8659 (3.7642 dilated) DO NOT USE THIS VALUE 

plot(time_surf(~isnan(float_sat_median_chl)),float_sat_median_chl(~isnan(float_sat_median_chl)),'*',time_surf(~isnan(float_sat_median_chl)),chl_surf(~isnan(float_sat_median_chl),1)/5.4315,'o');

