%{
 compare da output: old & new flux code versions
 8/20/15, BWB

 DYNAMO raw data computed with 2011 and 2015 versions of psd flux code
 Loads da files for each version and prepares plots

 2015 version da file format:
    COL     VARIABLE        DESCRIPTION
    1       doy             decimal day of year
                                                    2    3   4   5    6   7  8   9  10
    2:10    xship           COARE 3.5 input, scs:[Utrue,Tair,RH,P_mb,Tsea,Rs,Rl,lat,rain]
    11:15   no_data         To maintain historical column numbering...
                                                    16   17  18  19   20  21 22  23  24
    16:24   xpsd            COARE 3.5 input, psd:[Utrue,Tair,RH,P_mb,Tsea,Rs,Rl,lat,rain]
    25:29   no_data         To maintain historical column numbering...
    30      u_ship          ship mean true wind speed (m/s)
    31      u_ship_min      ship minimum true wind speed (m/s)
    32      u_ship_max      ship maximum true wind speed (m/s)
    33      u_ship_std      ship true wind speed stdev (m/s)
    34      wdir_ship       ship mean true wind direction
    35      wdir_ship_min   ship minimum true wind direction
    36      wdir_ship_max   ship maximum true wind direction
    37      wdir_ship_std   ship true wind direction stdev
    38      rwspd_psd       streamwise relative wind speed (m/s)from PSD sonic
    39      Un              N component, corrected true windspeed, PSD sonic
    40      Ue              E component, corrected true windspeed, PSD sonic
    41      wdir_psd        true wind direction from N,E components, PSD sonic
    42      rwdir_psd       relative wind direction from PSD sonic
    43      azm             azimuth angle from first rotation, radians
    44      tilt            tilt angle from first rotation, radians
    45      Tsea_psd        seasnake SST, degC
    46      Tair_psd        air temperature, PSD vaisala, degC
    47      RH_psd          RH, PSD vaisala
    48      RH_ship         RH, ship sensor
    49      Q_psd           specific humidity, PSD sensors, g/kg
    50      Q_ship          specific humidity, ship sensors, g/kg
    51      Qsea_psd        sea surface specific humidity from seasnake SST
    52      rhoa_psd         moist air density from PSD sensors (kg/m3)
                                                    53  54  55  56    57   58  59  60  61  62 63 64 65 66  67   68  69  70   71     72     73
    53:73   yship           COARE 3.5 output, ship:[usr,tau,hsb,hlb,hlwebb,tsr,qsr,zot,zoq,Cd,Ch,Ce,L,zet,dter,dqer,tkt,RF,Cdn_10,Chn_10,Cen_10]
                                                    74  75  76  77    78   79  80  81  82  83 84 85 86 87  88   89  90  91   92     93     94
    74:94   ypsd            COARE 3.5 output, PSD: [usr,tau,hsb,hlb,hlwebb,tsr,qsr,zot,zoq,Cd,Ch,Ce,L,zet,dter,dqer,tkt,RF,Cdn_10,Chn_10,Cen_10]
    95      Rl_std_psd      std PIR flux, PSD
    96      Rs_std_psd      std PSP flux, PSD
    97      sog_ship        mean sog, ship gps, m/s
    98      sog_ship_std    std sog, ship gps
    99      sog_ship_min    minimum sog, ship gps
    100     sog_ship_max    maximum sog, ship gps
    101     cog_ship        mean cog, deg, ship gps
    102     cog_ship_std    std cog, ship gps
    103     cog_ship_min    minimum cog, ship gps
    104     cog_ship_max    maximum cog, ship gps
    105     cog_psd         mean cog, deg, PSD gps
    106     cog_psd_std     std cog, PSD gps
    107     cog_psd_min     minimum cog, PSD gps
    108     cog_psd_max     maximum cog, PSD gps
    109     sog_psd         mean sog, m/s, PSD gps
    110     sog_psd_min     minimum sog, PSD gps
    111     sog_psd_max     maximum sog, PSD gps
    112     lrg             mean heading, deg, ship gyro
    113     lrg_std         std heading, ship gyro
    114     lrg_min         minimum heading, ship gyro
    115     lrg_max         maximum heading, ship gyro
    116     no_data         wu, Hare correction, not computed
    117     wu              wu covariance, Edson correction
    118     no_data         wv, Hare correction, not computed
    119     wv              wv covariance, Edson correction
    120     no_data         wT, Hare correction, not computed
    121     wt              wT covariance, Edson correction
    122     no_data         uvar, Hare correction, not computed
    123     no_data         vvar, Hare correction, not computed
    124     no_data         wvar, Hare correction, not computed
    125     uvar            u variance, Edson correction
    126     vvar            v variance, Edson correction
    127     wvar            w variance, Edson correction
    128     tvar            T variance
    129     Sp1             U platform velocity std, Edson correction
    130     Sp2             V platform velocity std, Edson correction
    131     Sp3             W platform velocity std, Edson correction
    132     hd_sin_std      stdev sine lrg heading
    133     hd_cos_std      stdev cosine lrg heading
    134     Cu2a            Cu structure func parameter, method a
    135     Cu2b            Cu structure func parameter, method b
    136     Cw2a            Cw structure func parameter, method a
    137     Cw2b            Cw structure func parameter, method b
    138     Ct2a            Ct structure func parameter, method a
    139     Ct2b            Ct structure func parameter, method b
    140     Cq2a            Ch2o structure func parameter, method a
    141     Cq2b            Ch2o structure func parameter, method b
    142     usib_a          ustar from idiss, method a
    143     usib_b          ustar from idiss, method b
    144     tsib_a          tstar from idiss, method a
    145     tsib_b          tstar from idiss, method b
    146     qsib_a          qstar from idiss, method a, kg/kg
    147     qsib_b          qstar from idiss, method b, kg/kg
    148     Lid_a           L from idiss, method a
    149     Lid_b           L from idiss, method b
    150     lon             longitude, scs gps
    151     lat             latitiude, scs gps
    152     Ubar            streamwise true windspeed from Edson correction
    153     odec            odec (ship speed through water), m/s
    154     odec_std        stdev odec
    155     odec_min        minimum odec
    156     odec_max        maximum odec
    157     P_mb            atmospheric pressure, PSD sensor, mb
extras:
    158     missing         missing data lines in son file per 10 min
    159     badPnts         # bad data points (spikes, etc.) in sonic data per 10 min
    160     hs_cov          sensible heat flux from covariance
    161     hl_cov          latent heat flux from covariance, uncorrected
    162     hlc_cov         webb corrected hl_cov
    163     sal_scs         salinity, ship system, ppt
    164     Ts_noise        Sonic temperature noise variance

 2011 differences: (no columns 160:164)
    COL     VARIABLE        DESCRIPTION         2     3    4   5  6  7  8  9   10  11  12  13 14 15
    2:15    xship           COARE 3.0 input, [U_true odec Tsea Ta Qs Qa Rs Rl rain zi P_mb zu zt zq]
                                                16    17   18  19 20 21 22 23  24  25  26  27 28 29
    16:29   xpsd            COARE 3.0 input, [U_true odec Tsea Ta Qs Qa Rs Rl rain zi P_mb zu zt zq]
                                              53  54  55  56 57  58  59 60 61  62  63   64   65  66   67    68 69 70   71     72     73
    53:73   yship           COARE 3.0 output,[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hl_webb Cd Ch Ce Cdn_10 Chn_10 Cen_10]
                                              74  75  76  77 78  79  80 81 82  83  84   85   86  87   88    89 90 91   92     93     94
    74:94   ypsd            COARE 3.0 output,[hsb hlb tau zo zot zoq L usr tsr qsr dter dqer tkt RF hl_webb Cd Ch Ce Cdn_10 Chn_10 Cen_10]
extras:
    158:159 slopeaccx       from motcorr routine??
%}

%% initialize
clear all
close all
savefigs = true; % set true to save png files
path_figs = '/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Compare_version_plots_c25';

%% load 2011 code output: N x 159 array
dfl = '/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_old/leg3/DYNAMO_2011_da_v3_vdecorr.txt';
v1 = load(dfl);

%% load 2015 code output: N x 164 array
dfl = '/Users/byronb/Raw_Data_Files.noindex/DYNAMO_2011/Revelle/flux/Processed_new/leg3/DYNAMO_2011_da_v3_decorr_c25.txt';
v3 = load(dfl);

%% bulk met plots

% true wind speed
figure; plot(v1(:,2),v3(:,2),'r.',v1(:,16),v3(:,16),'b.',[0,20],[0,20],'k'); grid;
axis([0,20,0,20]); title('True Wind Speed, v2015 vs v2011'); xlabel('U m/s, 2011 code');
ylabel('U m/s, 2015 code'); legend('SCS','PSD','location','northwest');

% true wind direction
figure; plot(v1(:,34),v3(:,34),'r.',v1(:,41),v3(:,41),'b.',[150,360],[150,360],'k'); grid;
axis([150,360,150,360]); title('True Wind Direction, v2015 vs v2011'); xlabel('dir, 2011 code');
ylabel('dir, 2015 code'); legend('SCS','PSD','location','northwest');

% Un
figure; plot(v1(:,39),v3(:,39),'b.',[-10,20],[-10,20],'k'); grid;
axis([-10,20,-10,20]); title('U north, v2015 vs v2011'); xlabel('U m/s, 2011 code');
ylabel('U m/s, 2015 code');

% Ue
figure; plot(v1(:,40),v3(:,40),'b.',[-5,20],[-5,20],'k'); grid;
axis([-5,20,-5,20]); title('U east, v2015 vs v2011'); xlabel('U m/s, 2011 code');
ylabel('U m/s, 2015 code');

% Tair
figure; plot(v1(:,5),v3(:,3),'r.',v1(:,19),v3(:,17),'b.',[0,35],[0,35],'k'); grid;
axis([23,30,23,30]); title('Air Temperature, v2015 vs v2011'); xlabel('Ta C, 2011 code');
ylabel('Ta C, 2015 code'); legend('SCS','PSD','location','northwest');

% Tsea
figure; plot(v1(:,4),v3(:,6),'r.',v1(:,18),v3(:,20),'b.',[0,35],[0,35],'k'); grid;
axis([28.5,30.5,28.5,30.5]); title('Sea Temperature, v2015 vs v2011'); xlabel('Ts C, 2011 code');
ylabel('Ts C, 2015 code'); legend('SCS','PSD','location','northwest');

% Rs
figure; plot(v1(:,8),v3(:,7),'r.',v1(:,22),v3(:,21),'b.',[0,1200],[0,1200],'k'); grid;
axis([0,1200,0,1200]); title('PSP, v2015 vs v2011'); xlabel('solar W/m2, 2011 code');
ylabel('solar W/m2, 2015 code'); legend('SCS','PSD','location','northwest');

% Rs std dev
figure; plot(v1(:,96),v3(:,96),'b.',[0,250],[0,250],'k'); grid;
axis([0,250,0,250]); title('PSD PSP std dev, v2015 vs v2011'); xlabel('solar std W/m2, 2011 code');
ylabel('solar std W/m2, 2015 code');

% Rl
figure; plot(v1(:,9),v3(:,8),'r.',v1(:,23),v3(:,22),'b.',[0,500],[0,500],'k'); grid;
axis([380,460,380,460]); title('PIR, v2015 vs v2011'); xlabel('IR W/m2, 2011 code');
ylabel('IR W/m2, 2015 code'); legend('SCS','PSD','location','northwest');

% Rl std dev
figure; plot(v1(:,95),v3(:,95),'b.',[0,12],[0,12],'k'); grid;
axis([0,12,0,12]); title('PSD PIR std dev, v2015 vs v2011'); xlabel('IR std W/m2, 2011 code');
ylabel('IR std W/m2, 2015 code');

% RH
figure; plot(v1(:,48),v3(:,48),'r.',v1(:,47),v3(:,47),'b.',[0,100],[0,100],'k'); grid;
axis([65,100,65,100]); title('Rel Humidity, v2015 vs v2011'); xlabel('RH, 2011 code');
ylabel('RH, 2015 code'); legend('SCS','PSD','location','northwest');

% Qa
figure; plot(v1(:,7),v3(:,50),'r.',v1(:,21),v3(:,49),'b.',[0,25],[0,25],'k'); grid;
axis([16,21,16,21]); title('Qair, v2015 vs v2011'); xlabel('Qa g/kg, 2011 code');
ylabel('Qa g/kg, 2015 code'); legend('SCS','PSD','location','northwest');

% Pressure
figure; plot(v1(:,157),v3(:,157),'b.',[1000,1015],[1000,1015],'k'); grid;
axis([1003,1010,1003,1010]); title('PSD Pressure, v2015 vs v2011'); xlabel('P mb, 2011 code');
ylabel('P mb, 2015 code');

% Relative wind spd
figure; plot(v1(:,38),v3(:,38),'b.',[0,20],[0,20],'k'); grid;
axis([0,20,0,20]); title('PSD Rel Wind Spd, v2015 vs v2011'); xlabel('Urel m/s, 2011 code');
ylabel('Urel m/s, 2015 code');

% Relative wind dir
ii = find(v1(:,42)>180); reldir = v1(:,42); reldir(ii) = reldir(ii)-360;
figure; plot(reldir,v3(:,42),'b.',[-90,90],[-90,90],'k'); grid;
axis([-90,90,-90,90]); title('PSD Rel Wind Dir, v2015 vs v2011'); xlabel('reldir, 2011 code');
ylabel('reldir, 2015 code');

% tilt
figure; plot(v1(:,44)*180/pi,v3(:,44)*180/pi,'b.',[0,10],[0,10],'k'); grid;
axis([3,7,3,7]); title('PSD tilt angle, v2015 vs v2011'); xlabel('tilt deg, 2011 code');
ylabel('tilt deg, 2015 code');

% gyro heading
figure; plot(v1(:,112),v3(:,112),'r.',[0,360],[0,360],'k'); grid;
axis([0,360,0,360]); title('Gyro Heading, v2015 vs v2011'); xlabel('lrg deg, 2011 code');
ylabel('lrg deg, 2015 code');

% rho air psd
figure; plot(v1(:,52),v3(:,52),'b.',[1,2],[1,2],'k'); grid;
axis([1.14,1.2,1.14,1.2]); title('PSD Air Density, v2015 vs v2011'); xlabel('Rho kg/m3, 2011 code');
ylabel('Rho kg/m3, 2015 code');

%wu
figure; plot(v1(:,117),v3(:,117),'b.',[-0.5,0.5],[-0.5,0.5],'k'); grid;
axis([-0.5,0.2,-0.5,0.2]); title('WU covariance, v2015 vs v2011'); xlabel('<wu>, 2011 code');
ylabel('<wu>, 2015 code');

%wv
figure; plot(v1(:,119),v3(:,119),'b.',[-0.5,0.5],[-0.5,0.5],'k'); grid;
axis([-0.2,0.2,-0.2,0.2]); title('WV covariance, v2015 vs v2011'); xlabel('<wv>, 2011 code');
ylabel('<wv>, 2015 code');

%wT
figure; plot(v1(:,121),v3(:,121),'b.',[-0.5,0.5],[-0.5,0.5],'k'); grid;
axis([-0.05,0.2,-0.05,0.2]); title('WTsonic covariance, v2015 vs v2011'); xlabel('<wT>, 2011 code');
ylabel('<wT>, 2015 code');

%uvar
figure; plot(v1(:,125),v3(:,125),'b.',[0,10],[0,10],'k'); grid;
axis([0,5,0,5]); title('U variance, v2015 vs v2011'); xlabel('uvar, 2011 code');
ylabel('uvar, 2015 code');

%vvar
figure; plot(v1(:,126),v3(:,126),'b.',[0,5],[0,5],'k'); grid;
axis([0,2,0,2]); title('V variance, v2015 vs v2011'); xlabel('vvar, 2011 code');
ylabel('vvar, 2015 code');

%wvar
figure; plot(v1(:,127),v3(:,127),'b.',[0,0.5],[0,0.5],'k'); grid;
axis([0,0.25,0,0.25]); title('W variance, v2015 vs v2011'); xlabel('wvar, 2011 code');
ylabel('wvar, 2015 code');

%tvar
figure; plot(v1(:,128),v3(:,128),'b.',[0,0.5],[0,0.5],'k'); grid;
axis([0,0.2,0,0.2]); title('Tsonic variance, v2015 vs v2011'); xlabel('tvar, 2011 code');
ylabel('tvar, 2015 code');

%sp1
figure; plot(v1(:,129),v3(:,129),'b.',[0,2],[0,2],'k'); grid;
axis([0,1,0,1]); title('Uplat std dev, v2015 vs v2011'); xlabel('std uplat, 2011 code');
ylabel('std uplat, 2015 code');

%sp2
figure; plot(v1(:,130),v3(:,130),'b.',[0,2],[0,2],'k'); grid;
axis([0,1,0,1]); title('Vplat std dev, v2015 vs v2011'); xlabel('std vplat, 2011 code');
ylabel('std vplat, 2015 code');

%sp3
figure; plot(v1(:,131),v3(:,131),'b.',[0,2],[0,2],'k'); grid;
axis([0,2,0,2]); title('Wplat std dev, v2015 vs v2011'); xlabel('std wplat, 2011 code');
ylabel('std wplat, 2015 code');

% ID Cu2
figure; plot(v1(:,134),v3(:,134),'r.',v1(:,135),v3(:,135),'b.',[0,1],[0,1],'k'); grid;
axis([0,0.25,0,0.25]); title('Inertial Dissipation Cu^2, v2015 vs v2011'); xlabel('Cu^2, 2011 code');
ylabel('Cu^2, 2015 code'); legend('Cu2a','Cu2b','location','northwest');

% ID Cw2
figure; plot(v1(:,136),v3(:,136),'r.',v1(:,137),v3(:,137),'b.',[0,1],[0,1],'k'); grid;
axis([0,0.25,0,0.25]); title('Inertial Dissipation Cw^2, v2015 vs v2011'); xlabel('Cw^2, 2011 code');
ylabel('Cw^2, 2015 code'); legend('Cw2a','Cw2b','location','northwest');

% ID Ct2
figure; plot(v1(:,138),v3(:,138),'r.',v1(:,139),v3(:,139),'b.',[0,1],[0,1],'k'); grid;
axis([0,0.02,0,0.02]); title('Inertial Dissipation Ct^2, v2015 vs v2011'); xlabel('Ct^2, 2011 code');
ylabel('Ct^2, 2015 code'); legend('Ct2a','Ct2b','location','northwest');

% ID Cq2
figure; plot(v1(:,140),v3(:,140),'r.',v1(:,141),v3(:,141),'b.',[0,1],[0,1],'k'); grid;
axis([0,0.05,0,0.05]); title('Inertial Dissipation Cq^2, v2015 vs v2011'); xlabel('Cq^2, 2011 code');
ylabel('Cq^2, 2015 code'); legend('Cq2a','Cq2b','location','northwest');

% ID usib
figure; plot(v1(:,142),v3(:,142),'r.',v1(:,143),v3(:,143),'b.',[0,1],[0,1],'k'); grid;
axis([0,0.8,0,0.8]); title('Inertial Dissipation usib, v2015 vs v2011'); xlabel('usib, 2011 code');
ylabel('usib, 2015 code'); legend('usib-a','usib-b','location','northwest');

% ID tsib
figure; plot(v1(:,144),v3(:,144),'r.',v1(:,145),v3(:,145),'b.',[-.2,.1],[-.2,.1],'k'); grid;
axis([-0.2,0.05,-0.2,0.05]); title('Inertial Dissipation tsib, v2015 vs v2011'); xlabel('tsib, 2011 code');
ylabel('tsib, 2015 code'); legend('tsib-a','tsib-b','location','northwest');

% ID qsib
figure; plot(v1(:,146),v3(:,146),'r.',v1(:,147),v3(:,147),'b.',[-1,1],[-1,1],'k'); grid;
axis([-0.003,0.001,-0.003,0.001]); title('Inertial Dissipation qsib, v2015 vs v2011'); xlabel('qsib, 2011 code');
ylabel('qsib, 2015 code'); legend('qsib-a','qsib-b','location','northwest');

% ID L
figure; plot(v1(:,148),v3(:,148),'r.',v1(:,149),v3(:,149),'b.',[-4000,0],[-4000,0],'k'); grid;
axis([-4000,0,-4000,0]); title('Inertial Dissipation L, v2015 vs v2011'); xlabel('lid m, 2011 code');
ylabel('lid m, 2015 code'); legend('lid-a','lid-b','location','northwest');

% u-star bulk
figure; plot(v1(:,60),v3(:,53),'r.',v1(:,81),v3(:,74),'b.',[-2,2],[-2,2],'k'); grid;
axis([0,0.7,0,0.7]); title('Ustar Bulk, v2015 vs v2011'); xlabel('usb m/s, 2011 C30 code');
ylabel('usb m/s, 2015 C35 code'); legend('ship','psd','location','northwest');

% stress bulk
figure; plot(v1(:,55),v3(:,54),'r.',v1(:,76),v3(:,75),'b.',[-2,2],[-2,2],'k'); grid;
axis([0,0.7,0,0.7]); title('Stress Bulk, v2015 vs v2011'); xlabel('tau, 2011 C30 code');
ylabel('tau, 2015 C35 code'); legend('ship','psd','location','northwest');

% t-star bulk
figure; plot(v1(:,61),v3(:,58),'r.',v1(:,82),v3(:,79),'b.',[-2,2],[-2,2],'k'); grid;
axis([-0.25,0.1,-0.25,0.1]); title('Tstar Bulk, v2015 vs v2011'); xlabel('tsb C, 2011 C30 code');
ylabel('tsb C, 2015 C35 code'); legend('ship','psd','location','northwest');

% q-star bulk
figure; plot(v1(:,62)*1000,v3(:,59),'r.',v1(:,83)*1000,v3(:,80),'b.',[-2,2],[-2,2],'k'); grid;
axis([-0.4,-0.1,-0.4,-0.1]); title('Qstar Bulk, v2015 vs v2011'); xlabel('qsb g/kg, 2011 C30 code');
ylabel('qsb g/kg, 2015 C35 code'); legend('ship','psd','location','northwest');

% hsb bulk
figure; plot(v1(:,53),v3(:,55),'r.',v1(:,74),v3(:,76),'b.',[-20,100],[-20,100],'k'); grid;
axis([0,100,0,100]); title('Hs Bulk, v2015 vs v2011'); xlabel('Hs W/m^2, 2011 C30 code');
ylabel('Hs W/m^2, 2015 C35 code'); legend('ship','psd','location','northwest');

% hlb bulk
figure; plot(v1(:,54),v3(:,56),'r.',v1(:,75),v3(:,77),'b.',[-20,500],[-20,500],'k'); grid;
axis([0,450,0,450]); title('Hl Bulk, v2015 vs v2011'); xlabel('Hl W/m^2, 2011 C30 code');
ylabel('Hl W/m^2, 2015 C35 code'); legend('ship','psd','location','northwest');

% Cool Skin delta-T
figure; plot(v1(:,63),v3(:,67),'r.',v1(:,84),v3(:,88),'b.',[-5,5],[-5,5],'k'); grid;
axis([0,0.4,0,0.4]); title('Cool Skin deltaT, v2015 vs v2011'); xlabel('dter C, 2011 C30 code');
ylabel('dter C, 2015 C35 code'); legend('ship','psd','location','northwest');


%% save all figures

if savefigs
    graphformat = '.png';
    graphdevice = '-dpng';
    figs = findall(0,'type','figure');
    sorted_figs = sortfighandlearray(figs); % put in ascending order
    N = length(sorted_figs);
    if N>0
        for ii=1:N
            figure(sorted_figs(ii));
            ppath = fullfile(path_figs,['fig_',sprintf('%02i',ii),graphformat]);
            print(graphdevice,ppath);
        end;
    end;
end;
