% plot PSD and SCS comparisons
% DEW 3/24/2008
% DEW only reads in SCS and ME2 data when ME1 was down

clear
close all
prtit = 1;

read_parameters_ICEALOT; %Reading parameter files flux conf from D:/data/lb/flux_eval_rhb_06
way_images_comp = 'S:\ICEALOT_2008\KNORR\Combined_instruments\comparisons\Raw_Images\';

% This is the cvs daily file and not the SCS from our system
md=input('Input mon and day to evaluate (0323 on).... ','s');
% This is the txt daily file and not the DAS raw file from our system
jd=input('Input yearday to evaluate (095 on).... ','s');

%Find all the Mean 1 txt files
% fn1=['*m1_' jd '.txt'];
% s1=[way_raw_data_flux fn1];
% dr1=dir(s1);
% [n1 m1]=size(dr1);
% mx1=n1;

way_raw_scs = 'S:\ICEALOT_2008\KNORR\SCS\Raw\';
%find all the scs txt files
fn2=['ICEALOT_scs_' md '.txt'];
s2=[way_raw_scs fn2];
%s1=['c:\data\icealot\data\Mailbox\' fname];
dr2=dir(s2);
[n2 m2]=size(dr2);
mx2=n2;

%Find all the Mean 2 txt files
fn3=['ICEALOT_m2_' jd '.txt'];
s3=[way_raw_data_flux fn3];
dr3=dir(s3);
[n3 m3]=size(dr3);
mx3=n3;


if mx1 ~= mx2
    display(['Number of daily text files not equal for mean and scs'])
    break
end
% DEW 3/27.2008  Put in to only do the last day
display('Only process latest day')
for i=mx2:mx2  %loop through temp profile text files
    % Needs to be 1 scs and 1 mean 1 and mean 2 for each day.
    clear m1
    %ICEALOT_M1_JJJ.txt;
    flnmmean1 = dr1(i).name
    f1=[way_raw_data_flux flnmmean1];
    ss1 = [f1];
    m1 = load(ss1);
    
    clear scs
    %ICEALOT_scs_MMDD.txt
    flnmscs = dr2(i).name
    f=[way_raw_scs flnmscs];
    ss2 = [f];
    scs = load(ss2);
    
    clear m2
    %ICEALOT_M1_JJJ.txt;
    flnmmean2 = dr3(i).name
    f3=[way_raw_data_flux flnmmean2];
    ss3 = [f3];
    m2 = load(ss3);
    
    
    if 1 % Mean 1 data
        org=m1(:,11); % optical rain gauge
        Tvais=m1(:,9);  %  Vaisala Temp
        Rhvais=m1(:,10);  % Vaisala RH
        Tsea=m1(:,8);   %correct for lead resistance 20C 2 wire
        %ii=find(Tvais<-50 || (Tvais>50 );Tvais(ii)=-999;    %ii=logical(Tvais<-100);Tvais(ii)=Nan;
    end
    
    %%%%%%%%%%%%
    if 1  % Mean 2 variables
        Tsea4 = m2(:,12);
        psp=m2(:,10);
        % % DEW Oct 14, 2006
        % PIR calculation
        sig_sb=5.67e-8;				        %Stefan Boltzmann constant
        Tc=m2(:,7);
        Td=m2(:,8);
        therm=m2(:,9);
        pir=therm+sig_sb*(Tc+273.16).^4-4*sig_sb*((Td+273.16).^4-(Tc+273.16).^4);
    end
    %%%%%%%%%%%%%%%
    jdd = str2num(jd);
    if jdd <= 80 % SCS data for 3/19-21 jd 79,80
        T1 = scs(:,1);
        T2 = scs(:,2);
        P1 = scs(:,3);
        cog = scs(:,7);
        sog = scs(:,8);
        P2 = scs(:,10);
        scssw = scs(:,11);
        rscs = scs(:,12);
        RH = scs(:,14);
        WD1 = scs(:,16);
        WS1 = scs(:,17);
        sst = scs(:,24);
        WDr = scs(:,26);
        WSr = scs(:,28);
        WDt = scs(:,27);
        WSt = scs(:,29);
    end
    
    if jdd <= 81 % SCS data for 3/19-3/21
%         T1 = scs(:,3);
%         T2 = scs(:,4);
%         P1 = scs(:,5);
%         cog = scs(:,9);
%         sog = scs(:,10);
%         P2 = scs(:,12);
%         scssw = scs(:,13);
%         rscs = scs(:,14);
%         RH = scs(:,16);
%         WD1 = scs(:,18);
%         WS1 = scs(:,19);
%         sst = scs(:,26);
%         WDr = scs(:,28);
%         WSr = scs(:,30);
%         WDt = scs(:,29);
%         WSt = scs(:,31);
        % From read_SCS_08a
         Ti = scs(:,2); %IMET_HRH AIR_TEMP
    Tv = scs(:,3); %WXT5_Ta Vaisala Weather Transmitter
    Pv = scs(:,4); %WXT5_Pa
    Lat = scs(:,5);
    Lon = scs(:,6);
    cog = scs(:,8);
    sog = scs(:,9);
    head = scs (:,10);
    Pi = scs(:,11);  % IMET
    scssw = scs(:,12);  % Eppley
    rscs = scs(:,13); %WXT rain intensity
    Rvc = scs(:,14); %WXT rain accumulation 
    RHi = scs(:,15); %IMET
    RHv = scs(:,16); % WXT RH
    WD1 = scs(:,17); %WXT Dir Relative
    WS1 = scs(:,18); %WXT Spd Relative
    SBsst = scs(:,21); % SeaBid SST
    TSGsst = scs(:,25); % TSG SST
    WDr = scs(:,27); % IMET RMY Dir True
    WDT = scs(:,28);  %WXT Dir True
    WSr = scs(:,29); % IMET RMY SPD True
    WST = scs(:,30); %WXT Spd True
       % 
    end
    if jdd >= 82  % daily files after 3/21
        Ti = scs(:,3); %IMET_HRH AIR_TEMP
        Tv = scs(:,4); %WXT5_Ta Vaisala Weather Transmitter
        Pv = scs(:,5); %WXT5_Pa
        Lat = scs(:,6);
        Lon = scs(:,7);
        cog = scs(:,9);
        sog = scs(:,10);
        head = scs (:,11);
        Pi = scs(:,12);  % IMET
        Esw = scs(:,13);  % Eppley
        Rvi = scs(:,14); %WXT rain intensity
        Rvc = scs(:,15); %WXT rain accumulation
        RHi = scs(:,16); %IMET RH
        RHv = scs(:,17); % WXT RH
        WDv = scs(:,18); %WXT Dir Relative
        WSv = scs(:,19); %WXT Spd Relative
        SBsst = scs(:,22); % SeaBid SST
        TSGsst = scs(:,26); % TSG SST
        WDrmy = scs(:,28); % IMET RMY Dir True
        WDvT = scs(:,29);  %WXT Dir True
        WSrmy = scs(:,30); % IMET RMY SPD True
        WSvT = scs(:,31); %WXT Spd True
    end
    [r,c] = size(m1);
    j1 = m1(:,c); %load julian day for mean 1 into array
    
    [r,c] = size(scs);
    j2 = scs(:,c); %load julian day from scs into array
    
    [r,c] = size(m2);
    j3 = m2(:,c); %load julian day from m2 into array
    
    jd = num2str(fix(j2(1,1)));
    if length(jd) == 2;
        jds(1) = '0';
        jds(2:3) = jd;
    end
    if length(jd) == 3
        jds(1:3) = jd;
    end
    s = way_images_comp;
    
    figure(1)
    set(gcf,'defaultaxesfonts',16)
    set(gca,'fontname','time','fontsize',16);
    plot(j2,Ti,'k',j2,Tv,'r',j1,Tvais,'b.');
    title('ICEALOT 2008');xlabel(['Year Day(UTC)']);ylabel('Temperature (C)');
    ylim([-15,15])
    legend('IMET','WXT','PSD','Location','Best')
    if prtit
        print_buffer = [s 'ICEALOT_PSDSCS_T_' jds '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(2)
    set(gcf,'defaultaxesfonts',16)
    set(gca,'fontname','time','fontsize',16);
    plot(j2,RHi,'k',j2,RHv,'r',j1,Rhvais,'b.');
    title('ICEALOT 2008');xlabel(['Year Day(UTC)']);ylabel('RH (%)');
    ylim([0,100])
    legend('IMET','WXT','PSD','Location','Best')
    if prtit
        print_buffer = [s 'ICEALOT_PSDSCS_RH_' jds '.png'];
        print('-dpng ', print_buffer);
    end
    
    figure(3)
    set(gcf,'defaultaxesfonts',16)
    set(gca,'fontname','time','fontsize',16);
    % plot(j2,sst,'k',j1,Tsea,'b.',j3,Tsea4,'r.');
    plot(j2,SBsst,'k',j2,TSGsst,'r',j3,Tsea4,'b.'); % skip 2 wire PSD probe  BROKEN
    title('ICEALOT 2008');xlabel(['Year Day(UTC)']);ylabel('SST (C)');
    msst = max(SBsst);;
    if msst > 20;
        msst = 15
    end
    ylim([-2,15])
    legend('Sea Bird','TSG','PSD','Location','Best')
    if prtit
        print_buffer = [s 'ICEALOT_PSDSCS_SST_' jds '.png'];
        print('-dpng ', print_buffer);
    end
    
    rsclr_08a  % calculate clear sky radiance
    
    figure(4);
    plot(j2,scssw,'k',j3,psp,'b.',j3,pir,'r.',j2,sw,'g');
    title('ICEALOT 2008');xlabel(['Year Day (UTC)']);ylabel('Solar FLux (W/m^2)');
    legend('SCS-PSP','PSD-PSP','PSD-PIR','Clear','Location','Best')
    ylim([-5 1000]);
    if prtit
        print_buffer = [s 'ICEALOT_PSDSCS_PSPPIR_' jds '.png'];
        print('-dpng ', print_buffer);
    end
    
end  % day loop