t% plot PSD and SCS comparisons
% DEW 3/24/2008

clear
close all

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


way_raw_mean = 'c:\data\icealot\data\mean\';
way_raw_scs = 'y:\knorr\scs\';

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

%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_mean fn3];
dr3=dir(s3);
[n3 m3]=size(dr3);
mx3=n3;

if mx2 ~= mx3 % SCS and mean 2
    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.
    if 0
        clear m1
        %ICEALOT_M1_JJJ.txt;
        flnmmean1 = dr1(i).name
        f1=[way_raw_mean 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);

    clear m2
    %ICEALOT_M2_JJJ.txt;
    flnmmean2 = dr3(i).name
    f3=[way_raw_mean flnmmean2];
    ss3 = [f3];
    m2 = load(ss3);

    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 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
    %%%%%%%%%%%%%%%
    if   0 % str2num(md(3:4)) < 22  %i < 4 % SCS data for 3/19-21
        str2num(md(3:4))
    Ti = scs(:,1); %IMET_HRH AIR_TEMP
    Tv = scs(:,2); %WXT5_Ta Vaisala Weather Transmitter
    Pv = scs(:,3); %WXT5_Pa
    Lat = scs(:,4);
    Lon = scs(:,5);
    cog = scs(:,7);
    sog = scs(:,8);
    head = scs (:,9);
    Pi = scs(:,10);  % IMET
    scssw = scs(:,11);  % Eppley
    Rvi = scs(:,12); %WXT rain intensity
    Rvc = scs(:,13); %WXT rain accumulation 
    RHi = scs(:,14); %IMET
    RHv = scs(:,15); % WXT RH
    WDv = scs(:,16); %WXT Dir Relative
    WSv = scs(:,17); %WXT Spd Relative
    SBsst = scs(:,20); % SeaBid SST
    TSGsst = scs(:,24); % TSG SST
    WDrmy = scs(:,26); % IMET RMY Dir True
    WDvT = scs(:,27);  %WXT Dir True
    WSrmy = scs(:,28); % IMET RMY SPD True
    WSvT = scs(:,29); %WXT Spd True
end

    if 1  %str2num(md(3:4)) >= 22  %i >= 4 % SCS data for 3/22 on
  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
    scssw = 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 fro mean 1 into array
    end

    [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(2:3) = jd;
            jds(1) = '0';  %make jd 3 characters
        end
        if length(jd) == 3
            jds = jd;
        end
    s = 'c:\data\icealot\data\comparisons\Raw_Images\'
if 0
    figure(1)
    set(gcf,'defaultaxesfonts',16)
    set(gca,'fontname','time','fontsize',16);
    plot(j2,T1,'k',j2,T2,'r',j3,Tsea4,'b.');
    title('ICEALOT 2008');xlabel(['Year Day(UTC)']);ylabel('Temperature (C)');
    ylim([-10,10])
    legend('SCS1','SCS2','PSD SST','Location','Best')
    print_buffer = [s 'ICEALOT_PSDSCS_T_' jds '.png'];
    print('-dpng ', print_buffer);

    figure(2)
        set(gcf,'defaultaxesfonts',16)
        set(gca,'fontname','time','fontsize',16);
        plot(j2,RH,'k',j1,Rhvais,'b.');
        title('ICEALOT 2008');xlabel(['Year Day(UTC)']);ylabel('RH (%)');
        ylim([0,100])
        legend('SCS','PSD','Location','Best')
        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,SBsst,'k',j2,TSGsst,'b',j3,Tsea4,'r.'); % 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('SCS Sea Bird','SCS TSG','PSD4','Location','Best')
    print_buffer = [s 'ICEALOT_PSDSCS_SST_' jds '.png'];

    print('-dpng ', print_buffer);


    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','Northwest')
    ylim([-5 1000]);
    print_buffer = [s 'ICEALOT_PSDSCS_PSPPIR_' jds '.png'];
    print('-dpng ', print_buffer);

    if 0
        figure(5)
        set(gcf,'defaultaxesfonts',16)
        set(gca,'fontname','time','fontsize',16);
        plot(j2,rscs,'k',j1,org,'r.');
        title('ICEALOT 2008');xlabel(['Year Day(UTC)']);ylabel('rain rate');
        ylim([-2,max(rscs)])
        legend('SCS','PSD','Location','Best')
        print_buffer = [s 'ICEALOT_PSDSCS_rain_' jds '.png'];
        print('-dpng ', print_buffer);
    end
end  % day loop