% plot only SCS comparisons of IMET and WXT when PSD was missing
% 90-102 (3/30-4/11)
% DEW 3/24/2008

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');

if 0  % skip psd data
    %Find all the Mean 1 txt files
    %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;
end

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;

if 0 % skip PSD data
    %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;
end

if 0  % do do test  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=1:1   %mx2:mx2  %loop through temp profile text files
    % Needs to be 1 scs and 1 mean 1 and mean 2 for each day.
    if 0
        clear m1
        %ICEALOT_M1_JJJ.txt;
        flnmmean1 = dr1(i).name
        f1=[way_raw_data_flux flnmmean1];
        ss1 = [f1];
        m1 = load(ss1);
    end
    
    clear scs
    %ICEALOT_scs_MMDD.txt
    flnmscs = dr2(i).name
    f=[way_raw_scs flnmscs];
    ss2 = [f];
    scs = load(ss2);
    if 0
        clear m2
        %ICEALOT_M1_JJJ.txt;
        flnmmean2 = dr3(i).name
        f3=[way_raw_data_flux flnmmean2];
        ss3 = [f3];
        m2 = load(ss3);
    end
    
    if 0 % 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 0  % 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);
        % From read_SCS_08a.m
    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 
    RH = 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 0  %jdd >= 82 % SCS data for 3/22 on
        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);
    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
    
    if 0
        [r,c] = size(m1);
        j1 = m1(:,c); %load julian day for mean 1 into array
        
        [r,c] = size(m2);
        j3 = m2(:,c); %load julian day from m2 into array
    end
    
    [r,c] = size(scs);
    j2 = scs(:,c); %load julian day from scs 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([-20,10])
    legend('IMET','WXT','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','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,msst])
    legend('Sea Bird','TSG','Location','Best')
    if prtit
        print_buffer = [s 'ICEALOT_PSDSCS_SST_' jds '.png'];
        print('-dpng ', print_buffer);
    end
    
    if 0 % skip radiometer plotting
        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
end  % day loop