%===============================================
%
% PLMR processing software
%
% This program calculates the brightness and physical temperatures
% from the PLMR counts. This software is meant to transfer counts
% to Tb and also as an internal (UoM/Chris) validation tool of 
% plmr2sp (ARA & Mahdi) and LMR2DAT (Arnstein), not as a stand alone program 
% for distribution!!! 
% This program is to be used after extracting the data from the cube with "cccp".
%
% Author          Date           version
% C. Rüdiger      02-Jan-09      1.0
% C. Rüdiger      04-Jan-09      1.01       added calibration run
% C. Rüdiger      18-Oct-09      2.0        converted into function
% C. Rüdiger      28-Nov-09      2.1        added calibration file reference
%
%===============================================

function f_raw_plmr

clear all
close all

%-----------------------------------
% import 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;

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

%-----------------------------------
% access calibration folder and source data files
%-----------------------------------
% load ProSensing factory calibrations
%cdpath=['cd ',fullfile(file_path,'calibrations')];
%eval(cdpath)
nm=['prosensing_',prosensing];
[a0v,a1v,a2v,a3v,a4v,a0h,a1h,a2h,a3h,a4h]=eval(genvarname(nm));

% load daily calibration files
if (calib==false)
    cdpath=['cd ',fullfile(file_path,file_date,'PLMRcal')];
    eval(cdpath)

    nm=['calib_',file_date(1:6),'_params'];
    [slpv,offsv,slph,offsh]=eval(genvarname(nm));
end

cdpath=['cd ',fullfile(file_path,file_date,'CUBE')];
eval(cdpath)

if (calib)
    file_name_calib=dir(fullfile(file_path,file_date,'PLMRcal',[file_date(1:6),'*.LMR']));
    sz=size(file_name_calib,1);
else
    sz=1;
end
    
tic     % start timer

for fi=1:sz
    
    if (calib)
        file_name_counts=file_name_calib(fi).name;
    end
    
    % counters
    i=0;
    c_T1=0;
    c_T2=0;
    c_T3=1;
    c_B=0;

    % initialise temperatures
    T3_calc(1,1:6)=40;
    tF=32;

    %-----------------------------------
    % open data files
    %-----------------------------------
    if (calib)
        fid=fopen(fullfile(file_path,file_date,'PLMRcal',file_name_counts));
    else
        fid=fopen(fullfile(file_path,file_date,'CUBE',file_name_counts));
    end

    %-----------------------------------
    % loop through the data file
    %-----------------------------------
    % pre-assign matrix (to save some computational time)
    Tbh=NaN(400000,7);
    Tbv=NaN(400000,7);
    T1_out=NaN(400000,7);
    T2_out=NaN(400000,7);
    T3_out=NaN(400000,7);

    %for i=1:25000
    while 1
        i=i+1;
        t = fgetl(fid);
        if ~ischar(t),   break,   end

        if (i==50000) disp('50000'); end
        if (i==100000) disp('100000'); end
        if (i==150000) disp('150000'); end
        if (i==200000) disp('200000'); end
        if (i==250000) disp('250000'); end

    %-----------------------------------
    % read in RADCON time stamp
    %-----------------------------------
        if (t(1,1)=='C')
            C=textscan(t,'%c %n %n *%n *%n *%n *%n');
            clk=(C{2}*65536+C{3})/512;

    %-----------------------------------
    % read in beam counts and transfer in gamma
    % and subsequently Tb (after the last beam);
    % dump current time stamp into the last column
    %-----------------------------------
    %    B=NaN(6,4);
    %    gamma_v=NaN(1,6);
    %    gamma_h=NaN(1,6);   
        elseif (t(1,1:3)=='B 0')
            c_B=c_B+1;
            in=str2double(t(1,3))+1;
            B(in,1)=str2double(t(5:9));
            B(in,2)=str2double(t(11:15));
            B(in,3)=str2double(t(17:21));
            B(in,4)=str2double(t(23:27));
            gamma_v(in)=(B(in,2)-B(in,3))/(B(in,1)-B(in,2));
            gamma_h(in)=(B(in,2)-B(in,4))/(B(in,1)-B(in,2));
        elseif (t(1,1:3)=='B 1')
            in=str2double(t(1,3))+1;
            B(in,1)=str2double(t(5:9));
            B(in,2)=str2double(t(11:15));
            B(in,3)=str2double(t(17:21));
            B(in,4)=str2double(t(23:27));
            gamma_v(in)=(B(in,2)-B(in,3))/(B(in,1)-B(in,2));
            gamma_h(in)=(B(in,2)-B(in,4))/(B(in,1)-B(in,2));
        elseif (t(1,1:3)=='B 2')
            in=str2double(t(1,3))+1;
            B(in,1)=str2double(t(5:9));
            B(in,2)=str2double(t(11:15));
            B(in,3)=str2double(t(17:21));
            B(in,4)=str2double(t(23:27));
            gamma_v(in)=(B(in,2)-B(in,3))/(B(in,1)-B(in,2));
            gamma_h(in)=(B(in,2)-B(in,4))/(B(in,1)-B(in,2));
        elseif (t(1,1:3)=='B 3')
            in=str2double(t(1,3))+1;
            B(in,1)=str2double(t(5:9));
            B(in,2)=str2double(t(11:15));
            B(in,3)=str2double(t(17:21));
            B(in,4)=str2double(t(23:27));
            gamma_v(in)=(B(in,2)-B(in,3))/(B(in,1)-B(in,2));
            gamma_h(in)=(B(in,2)-B(in,4))/(B(in,1)-B(in,2));
        elseif (t(1,1:3)=='B 4')
            in=str2double(t(1,3))+1;
            B(in,1)=str2double(t(5:9));
            B(in,2)=str2double(t(11:15));
            B(in,3)=str2double(t(17:21));
            B(in,4)=str2double(t(23:27));
            gamma_v(in)=(B(in,2)-B(in,3))/(B(in,1)-B(in,2));
            gamma_h(in)=(B(in,2)-B(in,4))/(B(in,1)-B(in,2));
        elseif (t(1,1:3)=='B 5')
            in=str2double(t(1,3))+1;
            B(in,1)=str2double(t(5:9));
            B(in,2)=str2double(t(11:15));
            B(in,3)=str2double(t(17:21));
            B(in,4)=str2double(t(23:27));
            gamma_v(1,in)=(B(in,2)-B(in,3))/(B(in,1)-B(in,2));
            gamma_h(1,in)=(B(in,2)-B(in,4))/(B(in,1)-B(in,2));
            Tbv(c_B,1:6)=a0v+a1v.*(tF/35)+a2v.*(T3_calc/35)+a3v.*gamma_v+a4v.*gamma_v.*(T3_calc/35);
            Tbh(c_B,1:6)=a0h+a1h.*(tF/35)+a2h.*(T3_calc/35)+a3h.*gamma_h+a4h.*gamma_h.*(T3_calc/35);
            if (calib==false)
                Tbv(c_B,1:6)=offsv'+slpv'.*Tbv(c_B,1:6);
                Tbh(c_B,1:6)=offsh'+slph'.*Tbh(c_B,1:6);
            end
            Tbv(c_B,7)=clk;
            Tbh(c_B,7)=clk;

    %-----------------------------------
    % read in temperature counts and
    % transfer into physical temperatures
    %-----------------------------------
        elseif (t(1,1:2)=='T1')
            T1(1)=str2double(t(4:6));
            T1(2)=str2double(t(8:10));
            T1(3)=str2double(t(12:14));
            T1(4)=str2double(t(16:18));
            T1(5)=str2double(t(20:22));
            T1(6)=str2double(t(24:26));
            c_T1=c_T1+1;
            T1_out(c_T1,1:6)=0.0548.*T1+6.60;
            T1_out(c_T1,7)=clk;
        elseif (t(1,1:2)=='T2')
            T2(1)=str2double(t(4:6));
            T2(2)=str2double(t(8:10));
            T2(3)=str2double(t(12:14));
            T2(4)=str2double(t(16:18));
            T2(5)=str2double(t(20:22));
            T2(6)=str2double(t(24:26));
            c_T2=c_T2+1;
            T2_out(c_T2,1:6)=0.0545.*T2+4.26;
            T2_out(c_T2,5)=0.0541.*T2(5)-6.89;       
            T2_out(c_T2,7)=clk;
            tF=mean(T2_out(c_T2,[1:4 6]));
        elseif (t(1,1:2)=='T3')
            T3(1)=str2double(t(4:8));
            T3(2)=str2double(t(10:14));
            T3(3)=str2double(t(16:20));
            T3(4)=str2double(t(22:26));
            T3(5)=str2double(t(28:32));
            T3(6)=str2double(t(34:38));
            T3_out(c_T3,1:6)=64.7-56.1.*(T3/65536)-20.9.*(T3/65536).^2;
            T3_calc=T3_out(c_T3,1:6);
            T3_out(c_T3,7)=clk;
            c_T3=c_T3+1;
        else
        end

    end    

    % close file and erase remaining unutilised matrix cells
    fclose(fid);
    Tbh(c_B+1:400000,:)=[];
    Tbv(c_B+1:400000,:)=[];
    T1_out(c_T1+1:400000,:)=[];
    T2_out(c_T2+1:400000,:)=[];
    T3_out(c_T3+1:400000,:)=[];

    %return

    %-----------------------------------
    % if this is a calibration run, calculate the mean for each polarization
    % and print the data into the corresponding data files
    %-----------------------------------
    if (calib)
        mean_Tbh=mean(Tbh(2:size(Tbh,1),1:6));
        mean_Tbv=mean(Tbv(2:size(Tbh,1),1:6));
        fidh=fopen(fullfile(file_path,file_date,'PLMRcal',['hpol_calib_',file_name_counts]),'w');
        fprintf(fidh,'%f,%f,%f,%f,%f,%f \n',Tbh(2:size(Tbh,1),1:6)');
        fidv=fopen(fullfile(file_path,file_date,'PLMRcal',['vpol_calib_',file_name_counts]),'w');
        fprintf(fidv,'%f,%f,%f,%f,%f,%f \n',Tbv(2:size(Tbv,1),1:6)');
        fclose('all');
    end

    %-----------------------------------
    % plot figures
    % Fig. 1: physical instrument temperatures
    % Fig. 2: brightness temperature time series
    % Fig. 3: brightness temperature scatter plots (3R vs 3L, 2R vs 2L, 1R vs 1L)
    %-----------------------------------
    figure(1)
    clf
    subplot(3,1,1)
    plot(T1_out(:,7),T1_out(:,1:5))
    ylabel('Temperature (Celsius)','fontsize',8);
    legend('Butler','Hot Load','Elec.BR','Elec.BL','Feed')
    subplot(3,1,2)
    plot(T2_out(:,7),T2_out(:,1:6))
    ylabel('Temperature (Celsius)','fontsize',8);
    legend('FL','BL','Mid','FR','Back','BR')
    subplot(3,1,3)
    plot(T3_out(:,7),T3_out(:,1:6))
    xlabel('Time','fontsize',8);ylabel('Temperature (Celsius)','fontsize',8);
    legend('3L','2L','1L','1R','2R','3R')
    if(plt)
        print('-djpeg','-r300',fullfile(file_path,file_date,'plots',['cal_physical_temperatures_',num2str(fi)]))
    end    

    figure(2)
    clf
    subplot(2,1,1)
    plot(Tbv(2:size(Tbh,1),7),Tbv(2:size(Tbh,1),1:6))
    title('Vertical Polarization')
    ylabel('Tb (Kelvin)','fontsize',8);
    legend('3L','2L','1L','1R','2R','3R')
    subplot(2,1,2)
    plot(Tbh(2:size(Tbh,1),7),Tbh(2:size(Tbh,1),1:6))
    title('Horizontal Polarization')
    legend('3L','2L','1L','1R','2R','3R')
    xlabel('Time','fontsize',8);ylabel('Tb (Kelvin)','fontsize',8);
    if(plt)
        print('-djpeg','-r300',fullfile(file_path,file_date,'plots',['cal_brightness_',num2str(fi)]))
    end

    figure(3)
    clf
    subplot(2,2,1)
    scatter(Tbv(2:size(Tbh,1),1),Tbv(2:size(Tbh,1),6),0.1,'r')
    hold on; box on
    scatter(Tbh(2:size(Tbh,1),1),Tbh(2:size(Tbh,1),6),0.1,'g')
    line([0 1000],[0 1000])
    xlim([100 320])
    ylim([100 320])
    xlabel('3L');ylabel('3R')

    subplot(2,2,2)
    scatter(Tbv(2:size(Tbh,1),2),Tbv(2:size(Tbh,1),5),0.1,'r')
    hold on; box on
    scatter(Tbh(2:size(Tbh,1),2),Tbh(2:size(Tbh,1),5),0.1,'g')
    line([0 1000],[0 1000])
    xlim([100 320])
    ylim([100 320])
    xlabel('2L');ylabel('2R')

    subplot(2,2,3)
    scatter(Tbv(2:size(Tbh,1),3),Tbv(2:size(Tbh,1),4),0.1,'r')
    hold on; box on
    scatter(Tbh(2:size(Tbh,1),3),Tbh(2:size(Tbh,1),4),0.1,'g')
    line([0 1000],[0 1000])
    xlim([100 320])
    ylim([100 320])
    xlabel('1L');ylabel('1R')
    if(plt)
        print('-djpeg','-r300',fullfile(file_path,file_date,'plots',['cal_brightness_temp_',num2str(fi)]))
    end

    fclose('all');

end

toc

end