%===============================================
%
% PLMR processing software (plotting)
%
% This program takes the data file from RASP
% and the RT3000 and creates various data plots. 
% This program is not a stand alone program for distribution!!! 
% It is to be used after processing the data with RASP.
%
% Author          Date           version
% C. R�diger      04-Jan-09      1.0
% C. R�diger      04-Jan-09      2.0     now reads in preprocessed
%                 data from RT and RASP output files
% C. R�diger      14-Sep-09      2.1     reads in new RASP format
% C. R�diger      12-Nov-09      3.0     changed to function
%
%===============================================

function f_graph_rasp

close all
clearvars %-except RASP_* var3_h

%-----------------------------------
% set parameters
%-----------------------------------

[linux,preproc,ussr,p2sp,rasp,calib,anci,kml,plmr_pro,graphs,experi,geb,pol,res,chk,sve,plt,tpause]=c_input_arg;
[fpath,file_date,file_name_PLMR,file_name_RASP,file_name_RT,file_name_counts, ... 
    prosensing,experi,pol,calib, ...
    plt] ... 
    =c_file_defs;
[high_res_start,high_res_end,low_res_start,low_res_end, ...
    minE,maxE,minN,maxN,UTM,datum,s_bin, ...
    temp_corr, ...
    RTstart,RTend,GPS_offset, ...
    step,range_h,bins_h,range_v,bins_v,range_TIR] ...
    =c_flight_params;

file_path=fullfile(fpath,'Data',experi,'field',file_date,'CUBE');

if (strcmp(res,'lo')==1)
    res2='low';
else
    res2='high';
end

%-----------------------------------
% open data files
%-----------------------------------
for j=1:size(pol,1)
    
    close all
    
    fname=fullfile(file_path,[file_name_RASP,'_',pol(j),'.dat']);
    fRASP=fopen(fname,'r');
    fgetl(fRASP);
    RASP=fscanf(fRASP,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
    RASP=reshape(RASP,15,size(RASP,1)/15)';
    fclose(fRASP);

%-----------------------------------
% temperature correction
%-----------------------------------
    RASP(:,6)=RASP(:,6)./(temp_corr(1).*(RASP(:,1)./86400./1000).^5 + ...
                          temp_corr(2).*(RASP(:,1)./86400./1000).^4 + ...
                          temp_corr(3).*(RASP(:,1)./86400./1000).^3 + ...
                          temp_corr(4).*(RASP(:,1)./86400./1000).^2 + ...
                          temp_corr(5).*(RASP(:,1)./86400./1000).^1 + ...
                          temp_corr(6));

%-----------------------------------
% save temperature corrected data
%-----------------------------------
    if(sve)
        fname=fullfile(file_path,['Rt_',file_name_RASP,'_',pol(j),'_',geb,'.dat']);
        fRASP=fopen(fname,'w');
        fprintf(fRASP,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',RASP');
        fclose(fRASP);
    end
                      
%-----------------------------------
% resample data according to beams and also geographic extend
% also filter all data from turns with more than a 2.5� roll
%-----------------------------------
    L3=RASP(RASP(:,5)==1,:);
    L2=RASP(RASP(:,5)==2,:);
    L1=RASP(RASP(:,5)==3,:);
    R1=RASP(RASP(:,5)==4,:);
    R2=RASP(RASP(:,5)==5,:);
    R3=RASP(RASP(:,5)==6,:);
    
    roll_f=2.5;

    if (strcmp(geb,'WH08')==1 || strcmp(geb,'WH09')==1)
        f_L3=find(L3(:,1)>low_res_start & L3(:,1)<low_res_end & ...
            L3(:,2)<550500 & abs(L3(:,12))<roll_f & L3(:,3)<maxN);
        f_L2=find(L2(:,1)>low_res_start & L2(:,1)<low_res_end & ...
            L2(:,2)<550500 & abs(L2(:,12))<roll_f & L2(:,3)<maxN);
        f_L1=find(L1(:,1)>low_res_start & L1(:,1)<low_res_end & ...
            L1(:,2)<550500 & abs(L1(:,12))<roll_f & L1(:,3)<maxN);
        f_R1=find(R1(:,1)>low_res_start & R1(:,1)<low_res_end & ...
            R1(:,2)<550500 & abs(R1(:,12))<roll_f & R1(:,3)<maxN);
        f_R2=find(R2(:,1)>low_res_start & R2(:,1)<low_res_end & ...
            R2(:,2)<550500 & abs(R2(:,12))<roll_f & R2(:,3)<maxN);
        f_R3=find(R3(:,1)>low_res_start & R3(:,1)<low_res_end & ...
            R3(:,2)<550500 & abs(R3(:,12))<roll_f & R3(:,3)<maxN);
    else
        f_L3=find(L3(:,1)>low_res_start & L3(:,1)<low_res_end & ... 
            abs(L3(:,12))<roll_f & L3(:,3)<maxN);
        f_L2=find(L2(:,1)>low_res_start & L2(:,1)<low_res_end & ...
            abs(L2(:,12))<roll_f & L2(:,3)<maxN);
        f_L1=find(L1(:,1)>low_res_start & L1(:,1)<low_res_end & ...
            abs(L1(:,12))<roll_f & L1(:,3)<maxN);
        f_R1=find(R1(:,1)>low_res_start & R1(:,1)<low_res_end & ...
            abs(R1(:,12))<roll_f & R1(:,3)<maxN);
        f_R2=find(R2(:,1)>low_res_start & R2(:,1)<low_res_end & ...
            abs(R2(:,12))<roll_f & R2(:,3)<maxN);
        f_R3=find(R3(:,1)>low_res_start & R3(:,1)<low_res_end & ...
            abs(R3(:,12))<roll_f & R3(:,3)<maxN);
    end

%-----------------------------------
% define starting points and pixel extend
%-----------------------------------
    minE=floor(minE/100)*100;
    maxE=ceil(maxE/100)*100;
    minN=floor(minN/100)*100;
    maxN=ceil(maxN/100)*100;

    E_start=minE+500;
    N_start=minN+500;

    E_diff=maxE-minE;
    N_diff=maxN-minN;

%-----------------------------------
% scale data for angular effects (reference currently beam 3L - outer left
% beam), assuming that they are correct and that the
% order of the data is the same throughout the file
% (ie. no data gaps or mix ups in the order of the beams)
%-----------------------------------
    %Tb_ref=mean([mean(var_s(1:6:size(var_s,1),3)) mean(var_s(6:6:size(var_s,1),3))]);
    Tb_ref=mean(L3(f_L3,6));

    L3_avg=mean(L3(f_L3,6));
    L2_avg=mean(L2(f_L2,6));
    L1_avg=mean(L1(f_L1,6));
    R1_avg=mean(R1(f_R1,6));
    R2_avg=mean(R2(f_R2,6));
    R3_avg=mean(R3(f_R3,6));

    L3(:,6)=L3(:,6).*(Tb_ref./L3_avg);
    L2(:,6)=L2(:,6).*(Tb_ref./L2_avg);
    L1(:,6)=L1(:,6).*(Tb_ref./L1_avg);
    R1(:,6)=R1(:,6).*(Tb_ref./R1_avg);
    R2(:,6)=R2(:,6).*(Tb_ref./R2_avg);
    R3(:,6)=R3(:,6).*(Tb_ref./R3_avg);

    L3_ref=L3(f_L3,:);
    L2_ref=L2(f_L2,:);
    L1_ref=L1(f_L1,:);
    R1_ref=R1(f_R1,:);
    R2_ref=R2(f_R2,:);
    R3_ref=R3(f_R3,:);

    y_lim(1)=floor((L1_avg-60)./10)*10;
    y_lim(2)=ceil((L1_avg+60)./10)*10;
%-----------------------------------
% split flight into north-south segments and remove variance
% by scaling to reference beam variance
%-----------------------------------

%-----------------------------------
% determine turns (only for Wirrangula Hill flights)
%-----------------------------------
    if (strcmp(geb,'WH08')==1 || strcmp(geb,'WH09')==1)
        fturn=0;
        for i=1:size(L3_ref,1)-1
            if (abs(L3_ref(i+1,13)-L3_ref(i,13))>45 & abs(L3_ref(i+1,13)-L3_ref(i,13))<270)
                %disp([num2str(L3_ref(i+1,13)),' ',num2str(L3_ref(i,13))]);
                fturn=fturn+1;
                tstamp(fturn)=i;
            end
        end

        sz_L3=1;
        sz_L2=1;
        sz_L1=1;
        sz_R1=1;
        sz_R2=1;
        sz_R3=1;

        clear Tb_L* Tb_R*

    %-----------------------------------
    % process individual transects
    %-----------------------------------
        for i=1:size(tstamp,2)+1
            if (i==1)
                filt_start=low_res_start;
                filt_end=L3_ref(tstamp(i),1);        
            elseif (1<i & i<=size(tstamp,2))
                filt_start=L3_ref(tstamp(i-1),1);
                filt_end=L3_ref(tstamp(i),1);
            else
                filt_start=filt_end;
                filt_end=low_res_end;
            end

            fL3=find(L3_ref(:,1)>filt_start & L3_ref(:,1)<filt_end);
            fL2=find(L2_ref(:,1)>filt_start & L2_ref(:,1)<filt_end);
            fL1=find(L1_ref(:,1)>filt_start & L1_ref(:,1)<filt_end);
            fR1=find(R1_ref(:,1)>filt_start & R1_ref(:,1)<filt_end);
            fR2=find(R2_ref(:,1)>filt_start & R2_ref(:,1)<filt_end);
            fR3=find(R3_ref(:,1)>filt_start & R3_ref(:,1)<filt_end);

            std_L3=std(L3_ref(fL3,6),1);
            std_L2=std(L2_ref(fL2,6),1);
            std_L1=std(L1_ref(fL1,6),1);
            std_R1=std(R1_ref(fR1,6),1);
            std_R2=std(R2_ref(fR2,6),1);
            std_R3=std(R3_ref(fR3,6),1);

            %m_L3=mean(L3_ref(fL3,6));
            %m_L2=mean(L2_ref(fL2,6));
            %m_L1=mean(L1_ref(fL1,6));
            %m_R1=mean(R1_ref(fR1,6));
            %m_R2=mean(R2_ref(fR2,6));
            %m_R3=mean(R3_ref(fR3,6));

            Tb_L3(sz_L3:sz_L3+size(fL3,1)-1,6)=L3_ref(fL3,6)+((1-(std_L3./std_L3))*(Tb_ref-L3_ref(fL3,6)));
            Tb_L2(sz_L2:sz_L2+size(fL2,1)-1,6)=L2_ref(fL2,6)+((1-(std_L3./std_L2))*(Tb_ref-L2_ref(fL2,6)));
            Tb_L1(sz_L1:sz_L1+size(fL1,1)-1,6)=L1_ref(fL1,6)+((1-(std_L3./std_L1))*(Tb_ref-L1_ref(fL1,6)));
            Tb_R1(sz_R1:sz_R1+size(fR1,1)-1,6)=R1_ref(fR1,6)+((1-(std_L3./std_R1))*(Tb_ref-R1_ref(fR1,6)));
            Tb_R2(sz_R2:sz_R2+size(fR2,1)-1,6)=R2_ref(fR2,6)+((1-(std_L3./std_R2))*(Tb_ref-R2_ref(fR2,6)));
            Tb_R3(sz_R3:sz_R3+size(fR3,1)-1,6)=R3_ref(fR3,6)+((1-(std_L3./std_R3))*(Tb_ref-R3_ref(fR3,6)));

            Tb_L3(sz_L3:sz_L3+size(fL3,1)-1,[1:5 7:15])=L3_ref(fL3,[1:5 7:15]);
            Tb_L2(sz_L2:sz_L2+size(fL2,1)-1,[1:5 7:15])=L2_ref(fL2,[1:5 7:15]);
            Tb_L1(sz_L1:sz_L1+size(fL1,1)-1,[1:5 7:15])=L1_ref(fL1,[1:5 7:15]);
            Tb_R1(sz_R1:sz_R1+size(fR1,1)-1,[1:5 7:15])=R1_ref(fR1,[1:5 7:15]);
            Tb_R2(sz_R2:sz_R2+size(fR2,1)-1,[1:5 7:15])=R2_ref(fR2,[1:5 7:15]);
            Tb_R3(sz_R3:sz_R3+size(fR3,1)-1,[1:5 7:15])=R3_ref(fR3,[1:5 7:15]);

            sz_L3=sz_L3+size(fL3,1);
            sz_L2=sz_L2+size(fL2,1);
            sz_L1=sz_L1+size(fL1,1);
            sz_R1=sz_R1+size(fR1,1);
            sz_R2=sz_R2+size(fR2,1);
            sz_R3=sz_R3+size(fR3,1);

        end

        L3=Tb_L3;
        L2=Tb_L2;
        L1=Tb_L1;
        R1=Tb_R1;
        R2=Tb_R2;
        R3=Tb_R3;
    
    end

    clear sz* Tb_L* Tb_R* fL* fR* std_* m_*

    var_2=[L3_ref;L2_ref;L1_ref;R1_ref;R2_ref;R3_ref];
    var_2=sortrows(var_2,1);

    if(sve)
%-----------------------------------
% save fully corrected data set to output file
%-----------------------------------
        fiz=fopen(fullfile(file_path,['R1_',file_name_RASP,'_',pol(j),'_',geb,'.dat']),'w');
        fprintf(fiz,'Lon, Lat, Tb, Time \n');
        fprintf(fiz,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',var_2');
        fclose('all')
    end

%-----------------------------------
% plot temperature-corrected Tb time series plots
%-----------------------------------
    figure(8)
    subplot(2,1,1)
    plot(RASP(RASP(:,5)==4,1)./1000,RASP(RASP(:,5)==4,6),'y')
    hold on
    plot(RASP(RASP(:,5)==5,1)./1000,RASP(RASP(:,5)==5,6),'r')
    plot(RASP(RASP(:,5)==6,1)./1000,RASP(RASP(:,5)==6,6),'k')
    title([geb,' ',pol(j),'-pol : Temperature Corrected T_b / Right hand beams']);
    xlabel('Seconds of Day since midnight UTC');
    ylabel('Corrected T_b [K]');ylim([200 300]);
    h=legend('R1','R2','R3');
    set(h);
    subplot(2,1,2)
    plot(RASP(RASP(:,5)==3,1)./1000,RASP(RASP(:,5)==3,6),'y')
    hold on
    plot(RASP(RASP(:,5)==2,1)./1000,RASP(RASP(:,5)==2,6),'r')
    plot(RASP(RASP(:,5)==1,1)./1000,RASP(RASP(:,5)==1,6),'k')
    title([geb,' ',pol(j),'-pol : Temperature Corrected T_b / Left hand beams']);
    xlabel('Seconds of Day since midnight UTC');
    ylabel('Corrected T_b [K]');ylim([200 300]);
    h=legend('L1','L2','L3');
    set(h);
    if(plt)
        print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['beams_raw_',pol(j),'_',geb]));
    end
    
%-----------------------------------
% plot temperature-corrected beam locations
%-----------------------------------
    figure(88)
    scatter(RASP(:,2),RASP(:,3),2,RASP(:,6),'filled')
    title(['Temperature corrected T_b [K] for ',geb,' in ',pol(j),'-pol'])
    xlabel('Easting [m]');ylabel('Northing [m]');box on;
    axis equal;box on;
    ylim([minN-6000 maxN+6000]);xlim([minE-6000 maxE+6000]);
    colorbar;caxis(eval(genvarname(['range_',pol(j)])));
    if(plt)
        print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['beams_raw_spatial_',pol(j),'_',geb]));
    end
    
%-----------------------------------
% plot fully corrected Tb time series plots
%-----------------------------------
    figure(9)
    subplot(2,1,1)
    plot(R1(:,1)./1000,R1(:,6),'y')
    hold on
    plot(R2(:,1)./1000,R2(:,6),'r')
    plot(R3(:,1)./1000,R3(:,6),'k')
    ylim(y_lim);
    title([geb,' ',pol(j),'-pol : Fully Corrected T_b / Right hand beams']);
    xlabel('Seconds of Day since midnight UTC');
    ylabel('Corrected T_b [K]');%ylim(eval(genvarname(['range_' pol(j)])))
    h=legend('R1','R2','R3');
    set(h);
    subplot(2,1,2)
    plot(L1(:,1)./1000,L1(:,6),'y')
    hold on
    plot(L2(:,1)./1000,L2(:,6),'r')
    plot(L3(:,1)./1000,L3(:,6),'k')
    ylim(y_lim);
    title([geb,' ',pol(j),'-pol : Fully Corrected T_b / Left hand beams']);
    xlabel('Seconds of Day since midnight UTC');
    ylabel('Corrected T_b [K]');%ylim(eval(genvarname(['range_' pol(j)])))
    h=legend('L1','L2','L3');
    set(h);
    if(plt)
        print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['beams_corr_',pol(j),'_',geb]))
    end

%-----------------------------------
% plot histograms of fully corrected Tb within previously defined range
%-----------------------------------
    figure(10)
    subplot(2,2,1)
    hist(var_2(:,6),eval(genvarname(['bins_' pol(j)])));
    xlim(eval(genvarname(['range_' pol(j)])));
    mu=mean(var_2(:,6));
    sig=std(var_2(:,6),1);
    %title(['Lowres all-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)])
    title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
    set(gca,'xticklabel',[]);
    subplot(2,2,2)
    hist([L3(:,6);R3(:,6)],eval(genvarname(['bins_' pol(j)])));
    xlim(eval(genvarname(['range_' pol(j)])));
    mu=mean([mean(L3(:,6));mean(R3(:,6))]);
    sig=std([L3(:,6);R3(:,6)],1);
    %title(['Lowres 3-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)])
    title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
    set(gca,'xticklabel',[]);
    subplot(2,2,3)
    hist([L2(:,6);R2(:,6)],eval(genvarname(['bins_' pol(j)])));    
    xlim(eval(genvarname(['range_' pol(j)])));
    mu=mean([mean(L2(:,6));mean(R2(:,6))]);
    sig=std([L2(:,6);R2(:,6)],1);
    %title(['Lowres 2-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)])
    title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
    xlabel('Brightness Temperature [K]');
    subplot(2,2,4)
    hist([L1(:,6);R1(:,6)],eval(genvarname(['bins_' pol(j)])));    
    xlim(eval(genvarname(['range_' pol(j)])));
    mu=mean([mean(L1(:,6));mean(R1(:,6))]);
    sig=std([L1(:,6);R1(:,6)],1);
    %title(['Lowres 1-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)])
    title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
    xlabel('Brightness Temperature [K]');
    if(plt)
        print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['lowres_hist_corr',pol(j),'_',geb]))
    end
    
%-----------------------------------
% plot histograms of fully corrected Tb for Lake Eyre and Simpson Desert
% (Nov 2008 data). This is necessary because of the distinct bi-modal data
% distribution at those two locations (and date).
%-----------------------------------
    if (strcmp(geb,'SD08')==1 || strcmp(geb,'LE08')==1 )
        if (strcmp(geb,'LE08')==1)
            sm1=45;
            lg1=130;
            sm2=220;
            lg2=311;
        else
            sm1=175;
            lg1=230;
            sm2=230;
            lg2=285;
        end        
        figure(100)
        subplot(2,2,[1 2])
        hist(var_2(:,6),[sm1:4:lg2])
        xlim([sm1+5 lg2-5]);
        mu=mean(var_2(:,6));
        sig=std(var_2(:,6),1);
    %    title(['Lowres all-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)])
        title(['Lowres all-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
        subplot(2,2,3)
        f=find(var_2(:,6)<lg1 & var_2(:,6)>sm1);
        hist(var_2(f,6),[sm1:4:lg1])
        xlim([sm1+5 lg1-5]);
        mu=mean(var_2(f,6));
        sig=std(var_2(f,6),1);
        title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
        subplot(2,2,4)
        f=find(var_2(:,6)>=sm2 & var_2(:,6)<lg2);
        hist(var_2(f,6),[sm2:4:lg2])
        xlim([sm2+5 lg2-5]);
        mu=mean(var_2(f,6));
        sig=std(var_2(f,6),1);
        title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8)
        if(plt)
            print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['lowres_hist2_',pol(j),'_',geb]))
        end
    end


%-----------------------------------
% start binning full data set onto a regular grid
%-----------------------------------

    Tb=NaN(N_diff/s_bin,E_diff/s_bin,3);

    clear *_avg

    for N=1:N_diff/s_bin
        if (N==25) disp('25'); end
        if (N==50) disp('50'); end
        if (N==75) disp('75'); end
        if (N==100) disp('100'); end
        if (N==250) disp('250'); end
        if (N==500) disp('500'); end
        if (N==750) disp('750'); end
        if (N==1000) disp('1000'); end

        for E=1:E_diff/s_bin

            E_center=E_start+(E-1)*s_bin;
            N_center=N_start+(N-1)*s_bin;

%-----------------------------------
% assume that all points within half the bin size from 
% the pixel's centre are contributing equally
% to the average and bin them into one pixel
%-----------------------------------
            num=find(var_2(:,2)>=E_center-s_bin/2 & var_2(:,2)<=E_center+s_bin/2 ...
                & var_2(:,3)>=N_center-s_bin/2 & var_2(:,3)<=N_center+s_bin/2);

            if (size(num,1)>0)
                Tb(N,E,1)=mean(var_2(num,6));
            end
            Tb(N,E,2)=E_center;
            Tb(N,E,3)=N_center;
        end    
    end

    r_Tb(:,1)=reshape(Tb(:,:,1),N_diff/s_bin*E_diff/s_bin,1);
    r_Tb(:,2)=reshape(Tb(:,:,2),N_diff/s_bin*E_diff/s_bin,1);
    r_Tb(:,3)=reshape(Tb(:,:,3),N_diff/s_bin*E_diff/s_bin,1);
    f=find(~isnan(r_Tb(:,1)));
    Tb_mean=mean(r_Tb(f,1));
    SD_Tb=std(r_Tb(f,1),1);

%-----------------------------------
% save gridded data to output file
%-----------------------------------
    if(sve)
        fiz=fopen(fullfile(file_path,['Tbgrid_',file_name_PLMR,'_',pol(j),'_',geb,'.dat']),'w');
        fprintf(fiz,'Tb,X_UTM,Y_UTM\n');
        fprintf(fiz,'%10.6f, %10.6f, %10.6f \n',r_Tb');
        fclose('all');
    end

%-----------------------------------
% plot spatial plot at 1km resolution
%-----------------------------------
    figure(1)
    clf
    pcolor(Tb(:,:,2),Tb(:,:,3),Tb(:,:,1))
    title([res2,' resolution Tb',pol(j),' [K] at ',geb])
    axis equal
    shading flat
    xlim([minE maxE]);%ylim([minN maxN]);
    colorbar; caxis(eval(genvarname(['range_' pol(j)])))
    xlabel('Easting [m]');ylabel('Northing [m]');
    if(plt)
        print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_Tb',pol(j),'_gridded_corr']))
    end

%-----------------------------------
% plot spatial anomalies plot at 1km resolution
%-----------------------------------
    figure(2)
    clf
    pcolor(Tb(:,:,2),Tb(:,:,3),Tb(:,:,1)-Tb_mean)
    title([res2,' resolution Tb',pol(j),' anomalies [K] at ',geb])
    axis equal
    shading flat
    xlim([minE maxE]);%ylim([minN maxN]);
    xlabel('Easting [m]');ylabel('Northing [m]');
    colorbar; caxis([-15*step 15*step])
    if(plt)
        print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_Tb',pol(j),'_gridded_anom']))
    end

%-----------------------------------
% plot histogram of 1km resolution data
%-----------------------------------
    figure(3)
    clf
    hist(r_Tb(f),[ceil(Tb_mean)-15*step:step:ceil(Tb_mean)+15*step])
    title([res2,' resolution Tb',pol(j),' at ',geb,' (mu=',num2str(Tb_mean), ...
        ', sigma=',num2str(SD_Tb),')'],'fontsize',12)
    %title(['mu=',num2str(Tb_mean),', sigma=',num2str(SD_Tb)],'fontsize',20,'color','r','fontweight','b')
    xlim([ceil(Tb_mean)-15*step ceil(Tb_mean)+15*step])
    set(gca,'fontsize',8)
    xlabel('Tb (Kelvin)');ylabel('# of occurences')
    if(plt)
        print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_Tb',pol(j),'_gridded_PDF']))
    end

    if (strcmp(geb,'SD08')==1 || strcmp(geb,'LE08')==1 || strcmp(geb,'LE09')==1)
        if (strcmp(geb,'LE08')==1)
            sm1=45;
            lg1=130;
            sm2=220;
            lg2=311;
        elseif (strcmp(geb,'SD08')==1)
            sm1=175;
            lg1=230;
            sm2=230;
            lg2=285;
        elseif (strcmp(geb,'LE09')==1)
            sm1=45;
            lg1=130;
            sm2=220;
            lg2=311;
        end        
        figure(103)
        subplot(2,2,[1 2])
        hist(r_Tb(f),[sm1:4:lg2])
        xlim([sm1+5 lg2-5]);
        mu=mean(r_Tb(f));
        sig=std(r_Tb(f),1);
        set(gca,'fontsize',8)
        ylabel('# of occurences','fontsize',8)
    %    title(['Lowres all-beams ',geb,', mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',20)
        title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',8,'color','r','fontweight','b')
        subplot(2,2,3)
        f=find(r_Tb<lg1 & r_Tb>sm1);
        hist(r_Tb(f),[sm1:4:lg1])
        xlim([sm1 lg1]);
        mu=mean(r_Tb(f));
        sig=std(r_Tb(f),1);
        set(gca,'xtick',[sm1:20:lg1],'fontsize',8)
        xlabel('Tb [Kelvin]','fontsize',8);ylabel('# of occurences','fontsize',8)
        title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',9,'color','k','fontweight','b')
        subplot(2,2,4)
        f=find(r_Tb>=sm2 & r_Tb<lg2);
        hist(r_Tb(f),[sm2:4:lg2])
        xlim([sm2 lg2]);
        mu=mean(r_Tb(f));
        sig=std(r_Tb(f),1);
        xlabel('Tb [Kelvin]','fontsize',8)
        set(gca,'xtick',[sm2:20:lg2],'fontsize',8)
        title(['mu=',num2str(mu),', sigma=',num2str(sig)],'fontsize',9,'color','k','fontweight','b')
        if(plt)
            print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',[geb,'_lowres_hist_gridded']))
        end
    end
end

fclose('all');

%-----------------------------------
% plot full scatter plot of low-resolution data
%-----------------------------------
fidh=fopen(fullfile(file_path,['Tbgrid_',file_name_PLMR,'_h_',geb,'.dat']),'r');
fidv=fopen(fullfile(file_path,['Tbgrid_',file_name_PLMR,'_v_',geb,'.dat']),'r');
fgetl(fidh);fgetl(fidv);
Tbgrid_h=fscanf(fidh,'%f,%f,%f\n',inf);
Tbgrid_h=reshape(Tbgrid_h,3,size(Tbgrid_h,1)./3);
Tbgrid_v=fscanf(fidv,'%f,%f,%f\n',inf);
Tbgrid_v=reshape(Tbgrid_v,3,size(Tbgrid_v,1)./3);
lat=Tbgrid_v(3,:);
lon=Tbgrid_v(2,:);

PR=(Tbgrid_v(1,:)-Tbgrid_h(1,:))./(Tbgrid_v(1,:)+Tbgrid_h(1,:));
v_h=Tbgrid_v(1,:)-Tbgrid_h(1,:);

figure(43)
subplot(1,2,1)
pcolor(reshape(lon,N_diff/s_bin,size(Tbgrid_v,2)/(N_diff/s_bin)) ...
    ,reshape(lat,N_diff/s_bin,size(Tbgrid_v,2)/(N_diff/s_bin)) ...
    ,reshape(PR,N_diff/s_bin,size(Tbgrid_v,2)/(N_diff/s_bin)))
axis equal;shading flat;box on;hold on;
title('Polarization Ratio [-]')
ylabel('Northing [m]');xlabel('Easting [m]');
clb=colorbar;set(clb,'Location','SouthOutside');
ylim([minN-1000 maxN+1000]);xlim([minE-1000 maxE+1000]);

subplot(1,2,2)
pcolor(reshape(lon,N_diff/s_bin,size(Tbgrid_v,2)/(N_diff/s_bin)) ...
    ,reshape(lat,N_diff/s_bin,size(Tbgrid_v,2)/(N_diff/s_bin)) ...
    ,reshape(v_h,N_diff/s_bin,size(Tbgrid_v,2)/(N_diff/s_bin)))
axis equal;shading flat;box on;hold on;
title('Polarization Difference [K]')
ylabel('Northing [m]');xlabel('Easting [m]');
clb=colorbar;set(clb,'Location','SouthOutside');
ylim([minN-1000 maxN+1000]);xlim([minE-1000 maxE+1000]);

if(plt)
    print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['PR_PD_',geb]));    
end

%-----------------------------------
% plot full scatter plot of low-resolution data
%-----------------------------------
fidh=fopen(fullfile(file_path,['Rt_',file_name_RASP,'_',pol(1),'_',geb,'.dat']),'r');
fidv=fopen(fullfile(file_path,['Rt_',file_name_RASP,'_',pol(2),'_',geb,'.dat']),'r');
fgetl(fidh);fgetl(fidv);

hpol=fscanf(fidh,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
hpol=reshape(hpol,15,size(hpol,1)/15)';
vpol=fscanf(fidv,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
vpol=reshape(vpol,15,size(vpol,1)/15)';

if (strcmp(geb,'WH08')==1 || strcmp(geb,'WH09')==1)
    fvpol=find(vpol(:,1)>low_res_start & vpol(:,1)<low_res_end & ...
        vpol(:,2)<550500 & abs(vpol(:,12))<roll_f & vpol(:,3)<maxN);
    fhpol=find(hpol(:,1)>low_res_start & hpol(:,1)<low_res_end & ...
        hpol(:,2)<550500 & abs(hpol(:,12))<roll_f & hpol(:,3)<maxN);
else
    fvpol=find(vpol(:,1)>low_res_start & vpol(:,1)<low_res_end & ... 
        abs(vpol(:,12))<roll_f & vpol(:,3)<maxN);
    fhpol=find(hpol(:,1)>low_res_start & hpol(:,1)<low_res_end & ... 
        abs(hpol(:,12))<roll_f & hpol(:,3)<maxN);
end

figure(4)
scatter(hpol(fhpol,7),hpol(fhpol,6),1,'filled','k');
hold on; box on;
scatter(vpol(fvpol,7),vpol(fvpol,6),1,'filled','r');
ylim(y_lim);
title(['Scatterplot of v & h pol data for ',geb]);
ylabel('T_b [K]');
xlabel('Incidence Angle [^o]');
h=legend('v-pol','h-pol');
set(h);
if(plt)
    print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['full_scatter_',geb]));    
end
fclose('all');

%-----------------------------------
% plot full scatter plot of low-resolution data (angular corrected)
%-----------------------------------
fidh=fopen(fullfile(file_path,['R1_',file_name_RASP,'_',pol(1),'_',geb,'.dat']),'r');
fidv=fopen(fullfile(file_path,['R1_',file_name_RASP,'_',pol(2),'_',geb,'.dat']),'r');
fgetl(fidh);fgetl(fidv);

hpol=fscanf(fidh,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
hpol=reshape(hpol,15,size(hpol,1)/15)';
vpol=fscanf(fidv,'%i,%f,%f,%f,%i,%f,%f,%i,%i,%i,%f,%f,%f,%f,%f\n',inf);
vpol=reshape(vpol,15,size(vpol,1)/15)';

if (strcmp(geb,'WH08')==1 || strcmp(geb,'WH09')==1)
    fvpol=find(vpol(:,1)>low_res_start & vpol(:,1)<low_res_end & ...
        vpol(:,2)<550500 & abs(vpol(:,12))<roll_f & vpol(:,3)<maxN);
    fhpol=find(hpol(:,1)>low_res_start & hpol(:,1)<low_res_end & ...
        hpol(:,2)<550500 & abs(hpol(:,12))<roll_f & hpol(:,3)<maxN);
else
    fvpol=find(vpol(:,1)>low_res_start & vpol(:,1)<low_res_end & ... 
        abs(vpol(:,12))<roll_f & vpol(:,3)<maxN);
    fhpol=find(hpol(:,1)>low_res_start & hpol(:,1)<low_res_end & ... 
        abs(hpol(:,12))<roll_f & hpol(:,3)<maxN);
end

figure(444)
scatter(hpol(fhpol,7),hpol(fhpol,6),1,'filled','k');
hold on; box on;
scatter(vpol(fvpol,7),vpol(fvpol,6),1,'filled','r');
ylim(y_lim);
title(['Scatterplot of v & h pol data for ',geb]);
ylabel('T_b [K]');
xlabel('Incidence Angle [^o]');
h=legend('v-pol','h-pol');
set(h);
if(plt)
    print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['full_scatter_corr_',geb]));    
end
fclose('all');

%-----------------------------------
% load RT data for plotting flight attitude information
%-----------------------------------
fid=fopen(fullfile(fpath,'Data',experi,'field',file_date,'RT3000', ...
    [file_name_RT,'.csv']),'r');
fgetl(fid);
RT=fscanf(fid,'%f/%f/%f,%f:%f:%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%i,\n',inf);
RT=reshape(RT,23,size(RT,1)/23)';
fclose('all');
RTtime=RT(:,7)+RTstart;

%-----------------------------------
% plot flight attitude information
%-----------------------------------
figure(44)
subplot(4,1,1)
plot(RTtime,RT(:,10))
ylabel('Altitude [m]');set(gca,'XTickLabel',[]);
xlim([min(RTtime) max(RTtime)]);
h=get(get(gca,'YLabel'),'Position');
title(['RT Flight Attitude Data for ',geb]);
subplot(4,1,2)
plot(RTtime,RT(:,18))
hold on;
line([0 1e6],[0 0],'color','k','linestyle','--')
ylabel('Pitch [^o]');set(gca,'XTickLabel',[]);
xlim([min(RTtime) max(RTtime)]);
h2=get(get(gca,'Ylabel'),'Position');
set(get(gca,'Ylabel'),'Position',[h(1);h2(2);h2(3)]);
ylim([-10 10])
subplot(4,1,3)
plot(RTtime,RT(:,19))
hold on;
line([0 1e6],[0 0],'color','k','linestyle','--')
ylabel('Roll [^o]');set(gca,'XTickLabel',[]);
xlim([min(RTtime) max(RTtime)]);
h2=get(get(gca,'Ylabel'),'Position');
set(get(gca,'Ylabel'),'Position',[h(1);h2(2);h2(3)]);
ylim([-20 20])
subplot(4,1,4)
plot(RTtime,RT(:,17))
ylabel('Heading [^o]');
h2=get(get(gca,'Ylabel'),'Position');
set(get(gca,'Ylabel'),'Position',[h(1);h2(2);h2(3)]);
xlabel('Time after GPS start [s]');
ylim([-10 370]);xlim([min(RTtime) max(RTtime)]);
if(plt)
    print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['flight_attitudes_',geb]));    
end

figure(4444)
plot(RTtime,RT(:,8))
ylabel('Latitude [^o]');
h2=get(get(gca,'Ylabel'),'Position');
set(get(gca,'Ylabel'),'Position',[h(1);h2(2);h2(3)]);
xlabel('Time after GPS start [s]');
ylim([-10 370]);xlim([min(RTtime) max(RTtime)]);


%-----------------------------------
% plot flight track
%-----------------------------------
figure(45)
subplot(2,1,1)
scatter(RT(:,9),RT(:,8),2,RT(:,10),'filled')
title(['Flight Track [^o] and Altitude ASL [m] for ',geb])
ylabel('Latitude [^o]');xlabel('Longitude [^o]');
t=colorbar;set(get(t,'title'),'String','Altitude [m]','HorizontalAlignment','left')
axis equal; box on;
subplot(2,1,2)
scatter(RT(:,9),RT(:,8),2,RT(:,17),'filled')
ylabel('Latitude [^o]');xlabel('Longitude [^o]');
t=colorbar;set(get(t,'title'),'String','Heading [^o]','HorizontalAlignment','left')
axis equal; box on;
if(plt)
    print('-depsc2','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['flight_track_',geb]));    
end


%-----------------------------------
% load electronics temperatures
%-----------------------------------
fid=fopen(fullfile(file_path,[file_name_PLMR,'_t.spa']),'r');
fgetl(fid);fgetl(fid);
T_out=fscanf(fid,'%i %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n',inf);
T_out=reshape(T_out,19,size(T_out,1)/19)';
f=find(T_out(:,1)>low_res_start & T_out(:,1)<low_res_end);
maxT1=ceil(max(max(T_out(:,2:6))))+1;
minT1=floor(min(min(T_out(:,2:6))))-1;
maxT2=ceil(max(max(T_out(:,8:13))))+1;
minT2=floor(min(min(T_out(:,8:13))))-1;
maxT3=ceil(max(max(T_out(:,14:19))))+1;
minT3=floor(min(min(T_out(:,14:19))))-1;

%-----------------------------------
% plot electronics temperatures --> only for low-res area
%-----------------------------------
figure(5)
set(gcf, 'PaperType','uslegal');
set(gcf, 'PaperUnits','normalized');
set(gcf, 'PaperOrientation','portrait');
%set(gcf, 'PaperSize', [1 2.5]);
subplot(3,1,1)
plot(T_out(f,1)./1000,T_out(f,2:6))
xlabel(' ','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
legend('Butler','Hot Load','Elec.BR','Elec.BL','Feed', ...
    'Location','SouthOutside','Orientation','horizontal')
ylim([32 42]);%([minT1 maxT1])
subplot(3,1,2)
plot(T_out(f,1)./1000,T_out(f,8:13))
xlabel(' ','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
legend('FL','BL','Mid','FR','Back','BR', ...
    'Location','SouthOutside','Orientation','horizontal')
ylim([20 35]);%ylim([minT2 maxT2])
subplot(3,1,3)
plot(T_out(f,1),T_out(f,14:19))
xlabel('Time','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
legend('3L','2L','1L','1R','2R','3R', ...
    'Location','SouthOutside','Orientation','horizontal')
ylim([35 45]);%ylim([minT3 maxT3])
if(plt)
    print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['physical_temperatures_',geb]))
end    

%-----------------------------------
% plot electronics temperatures --> all
%-----------------------------------
figure(6)
set(gcf, 'PaperType','uslegal');
set(gcf, 'PaperUnits','normalized');
set(gcf, 'PaperOrientation','portrait');
%set(gcf, 'PaperSize', [1 2.5]);
subplot(3,1,1)
plot(T_out(:,1)./1000,T_out(:,2:6))
xlabel(' ','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
legend('Butler','Hot Load','Elec.BR','Elec.BL','Feed', ...
    'Location','SouthOutside','Orientation','horizontal')
ylim([minT1 maxT1])
subplot(3,1,2)
plot(T_out(:,1)./1000,T_out(:,8:13))
xlabel(' ','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
legend('FL','BL','Mid','FR','Back','BR', ...
    'Location','SouthOutside','Orientation','horizontal')
ylim([minT2 maxT2])
subplot(3,1,3)
plot(T_out(:,1),T_out(:,14:19))
xlabel('Time','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
legend('3L','2L','1L','1R','2R','3R', ...
    'Location','SouthOutside','Orientation','horizontal')
ylim([minT3 maxT3])
if(plt)
    print('-djpeg','-r300',fullfile(fpath,'Data',experi,'field',file_date,'plots',['physical_temperatures_fulltime_',geb]))
end    


%===================================
%-----------------------------------
% quit script
%-----------------------------------
%===================================
return