disp('precip_comp3')
% read in raw SCS file and Mean1 txt file to compare OSI, IMET and WXT
% data only available until wave hit box. 079-089
clear
close all
fclose('all');
% This is the SCS raw daily file and not the SCS from our system
md=input('Input SCS raw mon and day to evaluate (0323).... ','s');
% This is the txt daily file and not the DAS raw file from our system
jd=input('Input MEAN1 yearday to evaluate (079-089 on).... ','s');

way_raw_mean = 'c:\data\icealot\data\mean\';
%Find Mean 2 txt files
fn2=['ICEALOT_m1_' jd '.txt'];
s3=[way_raw_mean fn2];
dr3=dir(s3);
[n3 m3]=size(dr3);
mx3=n3;
clear m1
%ICEALOT_M1_JJJ.txt;
flnmmean1 = dr3(1).name
f3=[way_raw_mean flnmmean1];
ss3 = [f3];
m1 = load(ss3);
[r,c] = size(m1);
org = m1(:,11); %PSD optical rain gauge
jdm1 = m1(:,c); % Year day PSD
jj = find(~isnan(org));
cumorg = cumsum(org(jj)); %cummulative precip
jdm = jdm1(jj);

ict = 0;
ict2 = 0;
way_raw_data = 'y:\Knorr\scs\';
e=[way_raw_data, 'KN08' md '.dat'];

disp(['Reading SCS data from raw file  ',e]);
fid=fopen(e,'r');

%R/V Knorr Calliope data file,  Sun 06/Apr/2008 00:00:00
%Current time zone: GMT Standard Time (GMT+0)
%Data timestamps are GMT
%****
%IMET_HRH  	39543.99996	23:59:56	77.669    0.254
%WXT5_Ta   	39543.99996	23:59:56	1
%WXT5_Ua   	39543.99996	23:59:56	75.6
%WXT5_Dm   	39543.99996	23:59:56	352
%WXT5_Sm   	39543.99996	23:59:56	17.6
%IMET_SWR  	39544.00000	23:59:59	0.4
%IMET_BPR  	39544.00000	23:59:59	1012.83
%PKEL99    	39544.00010	00:00:08	$PKEL99,000005,236.7,+005.00,00.00,+005.00,1500,70 32.051475N,018 26.944238E
%HEHDT     	39544.00011	00:00:09	$HEHDT,067.6,T*28
%GPS       	39544.00010	00:00:08	$GPZDA,000007.00,6,4,2008,+0,+0*69,$GPGGA,000008.00,7032.05131,N,01826.94953,E,2,8,1.4,26.94,M,32.37,M,10,0208*62,$GPVTG,115.2,T,,M,2.17,N,4.02,K*65
%IMET_PRC  	39544.00007	00:00:06	0.01    0.00  31.97
%TWind     	39544.00012	00:00:10	13.9276 67.5
dmp=fgetl(fid);		    %read 4 header lines
dmp=fgetl(fid);		    %read 4 header lines
dmp=fgetl(fid);		    %read 4 header lines
dmp=fgetl(fid);		    %read 4 header lines

%for p = 1:30  % short test read
while feof(fid)==0,
    %    clear st1 stz;
    %st1=fscanf(fid,'%g,%g,%g,%g,%g',[5,inf]);
    dtmp = fgetl(fid);
    if length(dtmp) > 33 & dtmp(1:8) == 'IMET_PRC'  %find IMET precip data
        ict = ict + 1;
        data1 = dtmp(33:length(dtmp));
        C = textscan(data1,'%f%f%f','CollectOutput',1);
        tmp = C{:};
        st1(2,ict) = tmp(1,1);  %
        st1(3,ict) = tmp(1,2);  %s
        st1(4,ict) = tmp(1,3);  %  rain total
        if 1
            %2008/03/22, 00:00:31
            tim = dtmp(24:31);
            YY = 2008;
            MM = str2num(md(:,1:2));
            DD = str2num(md(3:4));
            hh = str2num(tim(1:2));
            mm = str2num(tim(4:5));
            ss = str2num(tim(7:8));
            %calculate year day
            yd(ict,1) = datenum(YY,MM,DD,hh,mm,ss)-datenum(YY-1,12,31);
            year = num2str(YY);
        end
    end
    if length(dtmp) > 33 & dtmp(1:7) == 'WXT5_Rc'  %find IMET precip data
        ict2 = ict2 + 1;
        data1 = dtmp(33:length(dtmp));
        C = textscan(data1,'%f','CollectOutput',1);
        tmp = C{:};
        st2(2,ict2) = tmp(1,1);  %rain accumulation

        if 1
            %2008/03/22, 00:00:31
            tim = dtmp(24:31);
            YY = 2008;
            MM = str2num(md(:,1:2));
            DD = str2num(md(3:4));
            hh = str2num(tim(1:2));
            mm = str2num(tim(4:5));
            ss = str2num(tim(7:8));
            %calculate year day
            yd2(ict2,1) = datenum(YY,MM,DD,hh,mm,ss)-datenum(YY-1,12,31);
            year = num2str(YY);
        end
    end
end;                        %end while feof

if 1
    jdn = yd;
    jdn2 = yd2;
    s = 'c:\data\icealot\data\comparisons\Raw_Images\';
    figure(1);
    clf
    plot(jdn,st1(4,:),'k')
    hold
    plot(jdn2,st2(2,:),'r');
    title('Knorr SCS ')
    xlabel('Year Day (UTC)');
    ylabel('rain (mm)');
    mm(1) = max(st1(4,:));
    mm(2) = max(st2(2,:));
    mr = max(mm);
    ylim([0 mr+5])
    legend('IMET', 'WXT','Location','Best')
    %print_buffer = [s 'ICEALOT_SCS_rain_' md '.png'];
    %print('-dpng ', print_buffer);

    figure(2);
    clf
    plot(jdm,cumorg/10);
    title('PSD org rain rate');
    xlabel('Year Day (UTC)');
    ylabel('Rain(cm)');

    %subtract the minimum from the Imet sensor to bring it into scale for plotting
    IMETr = st2(2,:) - min(st2(2,:));
    figure(3);
    clf
    plot(jdn,st1(4,:),'k.')
    hold
    plot(jdn2,IMETr,'r.');
    plot(jdm,cumorg/10,);
    title('ICEALOT 2008 ')
    xlabel('Year Day (UTC)');
    ylabel('Rain');
    mm(1) = max(st1(4,:));
    mm(2) = max(IMETr)+.5;
    mr = max(mm);
    %ylim([0 mr+5])
    legend('IMET', 'WXT','PSD-OSI','Location','Best')
    print_buffer = [s 'ICEALOT_rain_' md '.png'];
    print('-dpng ', print_buffer);
end %skip any plots