%===============================================
%
% PLMR processing software (calibration)
%
% This program takes the data files from plmr.m
% and calculates the slope and offset for the calibration
% parameters. This is not a stand alone program for 
% distribution!!! This program is to be used after
% processing the data with plmr.m.
%
% Don't forget to change the comment lines in the output files!
%
% Author          Date           version
% C. Rüdiger      04-Jan-09      1.0
% C. Rüdiger      17-Jan-09      1.1    added flight days and statistics
% C. Rüdiger      18-Oct-09      2.0    converted into function
% C. Rüdiger      28-Nov-09      2.1    changed to calibration to search
%                                       for Tb files created with f_raw_plmr
% C. Rüdiger      18-Mar-10      2.11   corrected wrong file reference
%
%===============================================

function f_calibration

clear all
close all

%-----------------------------------
% 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;
calib_path=fullfile(fpath,'Data',experi,'field',file_date,'PLMRcal');

sky_temp=8.;

fid=dir(fullfile(calib_path,[file_date(1:6),'*Pre.txt']));
pre_temp=fid.name;
fid=dir(fullfile(calib_path,[file_date(1:6),'*Post.txt']));
post_temp=fid.name;
fid=dir(fullfile(calib_path,['*pol*',file_date(1:6),'*PreWarm.LMR']));
pre_warm_h=fid(1).name;
pre_warm_v=fid(2).name;
fid=dir(fullfile(calib_path,['*pol*',file_date(1:6),'*PostWarm.LMR']));
post_warm_h=fid(1).name;
post_warm_v=fid(2).name;
fid=dir(fullfile(calib_path,['*pol*',file_date(1:6),'*PreCold.LMR']));
pre_cold_h=fid(1).name;
pre_cold_v=fid(2).name;
fid=dir(fullfile(calib_path,['*pol*',file_date(1:6),'*PostCold.LMR']));
post_cold_h=fid(1).name;
post_cold_v=fid(2).name;

%-----------------------------------
% load data (vpol & hpol)
%-----------------------------------
pre_warm_vpol=load(fullfile(calib_path,[pre_warm_v]));
pre_warm_hpol=load(fullfile(calib_path,[pre_warm_h]));
pre_cold_vpol=load(fullfile(calib_path,[pre_cold_v]));
pre_cold_hpol=load(fullfile(calib_path,[pre_cold_h]));
pre_temp_box=load(fullfile(calib_path,[pre_temp]));

post_warm_vpol=load(fullfile(calib_path,[post_warm_v]));
post_warm_hpol=load(fullfile(calib_path,[post_warm_h]));
post_cold_vpol=load(fullfile(calib_path,[post_cold_v]));
post_cold_hpol=load(fullfile(calib_path,[post_cold_h]));
post_temp_box=load(fullfile(calib_path,[post_temp]));

%-----------------------------------
% prepare data (vpol & hpol)
%-----------------------------------
pre_box=mean(mean(pre_temp_box))+273.15;
post_box=mean(mean(post_temp_box))+273.15;

mean_pre_warm_vpol=mean(pre_warm_vpol);
mean_pre_warm_hpol=mean(pre_warm_hpol);
mean_pre_cold_vpol=mean(pre_cold_vpol);
mean_pre_cold_hpol=mean(pre_cold_hpol);

mean_post_warm_vpol=mean(post_warm_vpol);
mean_post_warm_hpol=mean(post_warm_hpol);
mean_post_cold_vpol=mean(post_cold_vpol);
mean_post_cold_hpol=mean(post_cold_hpol);

%-----------------------------------
% compute statistics for the observed
% brightness temperatures (vpol & hpol)
%-----------------------------------
SD_pre_warm_vpol=std(pre_warm_vpol,1);
SD_pre_warm_hpol=std(pre_warm_hpol,1);
SD_pre_cold_vpol=std(pre_cold_vpol,1);
SD_pre_cold_hpol=std(pre_cold_hpol,1);

SD_post_warm_vpol=std(post_warm_vpol,1);
SD_post_warm_hpol=std(post_warm_hpol,1);
SD_post_cold_vpol=std(post_cold_vpol,1);
SD_post_cold_hpol=std(post_cold_hpol,1);

%-----------------------------------
% compute slope and offset (vpol & hpol)
%-----------------------------------
% pre-flight calibration
%-----------------------------------
slp_vpol_pre=(pre_box-sky_temp)./(mean_pre_warm_vpol-mean_pre_cold_vpol);
off_vpol_pre=sky_temp-slp_vpol_pre.*mean_pre_cold_vpol;
slp_hpol_pre=(pre_box-sky_temp)./(mean_pre_warm_hpol-mean_pre_cold_hpol);
off_hpol_pre=sky_temp-slp_hpol_pre.*mean_pre_cold_hpol;

%-----------------------------------
% post-flight calibration
%-----------------------------------
slp_vpol_post=(post_box-sky_temp)./(mean_post_warm_vpol-mean_post_cold_vpol);
off_vpol_post=sky_temp-slp_vpol_post.*mean_post_cold_vpol;
slp_hpol_post=(post_box-sky_temp)./(mean_post_warm_hpol-mean_post_cold_hpol);
off_hpol_post=sky_temp-slp_hpol_post.*mean_post_cold_hpol;

slp_vpol=(slp_vpol_pre+slp_vpol_post)./2.
slp_hpol=(slp_hpol_pre+slp_hpol_post)./2.
off_vpol=(off_vpol_pre+off_vpol_post)./2.
off_hpol=(off_hpol_pre+off_hpol_post)./2.

%-----------------------------------
% verify calibration by calculating target temperatures 
% (warm and cold references)
%-----------------------------------
target_pre_warm=slp_vpol_pre.*mean_pre_warm_vpol+off_vpol_pre
target_pre_cold=slp_vpol_pre.*mean_pre_cold_vpol+off_vpol_pre
target_post_warm=slp_vpol_post.*mean_post_warm_vpol+off_vpol_post
target_post_cold=slp_vpol_post.*mean_post_cold_vpol+off_vpol_post

scat_pre_box=pre_box*ones(1,6);
scat_post_box=post_box*ones(1,6);
scat_sky_temp=sky_temp*ones(1,6);

%-----------------------------------
% create output figures
%-----------------------------------
figure(1)
clf
subplot(2,2,1)
scatter(mean_pre_warm_vpol,scat_pre_box,1.5,'filled')
title('Pre Vpol')
hold on; box on;
scatter(mean_pre_cold_vpol,scat_sky_temp,1.5,'filled')
plot([mean_pre_cold_vpol;mean_pre_warm_vpol],[scat_sky_temp;scat_pre_box])
xlabel('PLMR');ylabel('Obs');
xlim([-20 320]);ylim([0 320]);

subplot(2,2,2)
scatter(mean_pre_warm_hpol,scat_pre_box,1.5,'filled')
title('Pre Hpol')
hold on; box on;
scatter(mean_pre_cold_hpol,scat_sky_temp,1.5,'filled')
plot([mean_pre_cold_hpol;mean_pre_warm_hpol],[scat_sky_temp;scat_pre_box])
xlabel('PLMR');ylabel('Obs');
xlim([-20 320]);ylim([0 320]);

subplot(2,2,3)
scatter(mean_post_warm_vpol,scat_post_box,1.5,'filled')
title('Post Vpol')
hold on; box on;
scatter(mean_post_cold_vpol,scat_sky_temp,1.5,'filled')
plot([mean_post_cold_vpol;mean_post_warm_vpol],[scat_sky_temp;scat_post_box])
xlabel('PLMR');ylabel('Obs');
xlim([-20 320]);ylim([0 320]);

subplot(2,2,4)
scatter(mean_post_warm_hpol,scat_post_box,1.5,'filled')
title('Post Hpol')
hold on; box on;
scatter(mean_post_cold_hpol,scat_sky_temp,1.5,'filled')
plot([mean_post_cold_hpol;mean_post_warm_hpol],[scat_sky_temp;scat_post_box])
xlabel('PLMR');ylabel('Obs');
xlim([-20 320]);ylim([0 320]);

if (plt)
    print('-djpeg','-r300',fullfile(calib_path,[file_date,'_calib']));
end

figure(2)
clf
subplot(2,2,1)
plot(pre_warm_vpol)
title('Pre Vpol warm')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);

subplot(2,2,2)
plot(pre_cold_vpol)
title('Pre Vpol cold')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);
legend('3L','2L','1L','1R','2R','3R')

subplot(2,2,3)
plot(post_warm_vpol)
title('Post Vpol warm')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);

subplot(2,2,4)
plot(post_cold_vpol)
title('Post Vpol cold')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);

if (plt)
    print('-djpeg','-r300',fullfile(calib_path,[file_date,'_calib_v']));
end

figure(3)
clf
subplot(2,2,1)
plot(pre_warm_hpol)
title('Pre Hpol warm')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);

subplot(2,2,2)
plot(pre_cold_hpol)
title('Pre Hpol cold')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);
legend('3L','2L','1L','1R','2R','3R')

subplot(2,2,3)
plot(post_warm_hpol)
title('Post Hpol warm')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);

subplot(2,2,4)
plot(post_cold_hpol)
title('Post Hpol cold')
hold on; box on;
xlabel('time');ylabel('Obs');
%xlim([-20 320]);ylim([0 320]);

if (plt)
    print('-djpeg','-r300',fullfile(calib_path,[file_date,'_calib_h']));
end

%-----------------------------------
% create output files
%-----------------------------------
if (chk)
    fid=fopen(fullfile(calib_path,[file_date,'_params_chk.m']),'w');
    fprintf(fid,['function [slpvpre,offsvpre,slphpre,offshpre,slpvpost,' ...
        'offsvpost,slphpost,offshpost]=calib_params\n\n']);
    fprintf(fid,'slpvpre=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',slp_vpol_pre);
    fprintf(fid,'offsvpre=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',off_vpol_pre);
    fprintf(fid,'slphpre=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',slp_hpol_pre);
    fprintf(fid,'offshpre=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n\n',off_hpol_pre);
    fprintf(fid,'slpvpost=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',slp_vpol_post);
    fprintf(fid,'offsvpost=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',off_vpol_post);
    fprintf(fid,'slphpost=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',slp_hpol_post);
    fprintf(fid,'offshpost=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n\n',off_hpol_post);    
    fprintf(fid,'end');
    fclose(fid);
    
    fid=fopen(fullfile(calib_path,['calib_',file_date,'_params.m']),'w');
    fprintf(fid,'function [slpv,offsv,slph,offsh]=calib_params\n\n');
    fprintf(fid,'slpv=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',slp_vpol);
    fprintf(fid,'offsv=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',off_vpol);
    fprintf(fid,'slph=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n',slp_hpol);
    fprintf(fid,'offsh=[%8.6f;%8.6f;%8.6f;%8.6f;%8.6f;%8.6f];\n\n',off_hpol);
    fprintf(fid,'end');
    fclose(fid);

    fid=fopen(fullfile(calib_path,[file_date,'_plmr.cal']),'w');
    fprintf(fid,'[calibration]\n');
    
    if (strcmp(geb,'WH09')==1 || strcmp(geb,'SD09')==1 || strcmp(geb,' P01')==1 || ...
            strcmp(geb,' P02')==1 || strcmp(geb,' P03')==1 || strcmp(geb,' P04')==1 || ... 
            strcmp(geb,' P05')==1 || strcmp(geb,' P06')==1 || strcmp(geb,' P07')==1 || ...
            strcmp(geb,' P08')==1 || strcmp(geb,' P09')==1 || strcmp(geb,' P10')==1 || ...
            strcmp(geb,'1T00')==1 || strcmp(geb,'2T00')==1 || strcmp(geb,'3T00')==1)
        fprintf(fid,'Date=2009.07.17\n');
        fprintf(fid,'Description=pre-SAZEOz-II (July 2009 ff.)\n');
        fprintf(fid,'Comment=ProSensing calibration\n');
        fprintf(fid,'Parameters=Offset,Antenna_Mean,Receiver,Gamma,Gamma*Receiver\n');
        fprintf(fid,'Beam3LV= 295.89  -26.09   48.51  256.67  -16.70\n');
        fprintf(fid,'Beam2LV= 283.51  -22.62   66.62  236.51    5.92\n');
        fprintf(fid,'Beam1LV= 261.21  -16.61   70.02  217.98   33.48\n');
        fprintf(fid,'Beam1RV= 268.27  -18.53   63.37  226.02  -19.16\n');
        fprintf(fid,'Beam2RV= 290.97  -16.82   43.51  239.18  -28.70\n');
        fprintf(fid,'Beam3RV= 281.81  -12.50   45.40  225.26  -16.60\n');
        fprintf(fid,'Beam3LH= 296.47  -33.37   53.51  260.52  -25.84\n');
        fprintf(fid,'Beam2LH= 279.68  -31.86   74.86  228.22    2.66\n');
        fprintf(fid,'Beam1LH= 269.96  -29.93   79.78  216.12   33.43\n');
        fprintf(fid,'Beam1RH= 273.69  -30.02   69.20  224.26  -19.76\n');
        fprintf(fid,'Beam2RH= 289.52  -26.91   54.17  233.09  -26.18\n');
        fprintf(fid,'Beam3RH= 282.64  -31.13   61.62  221.48  -14.97\n\n');
    else
        fprintf(fid,'Date=2008.09.09\n');
        fprintf(fid,'Description=pre-Tipperary (2008)\n');
        fprintf(fid,'Comment=ProSensing calibration using data after fixing PLMR\n');
        fprintf(fid,'Parameters=Offset,Antenna_Mean,Receiver,Gamma,Gamma*Receiver\n');
        fprintf(fid,'Beam3LV= 307.53  -13.75   29.67  179.23   63.91\n');
        fprintf(fid,'Beam2LV= 275.82  -12.87   57.65  136.04   83.05\n');
        fprintf(fid,'Beam1LV= 278.66  -13.97   45.27  229.05   20.01\n');
        fprintf(fid,'Beam1RV= 243.33  -76.41  133.81  101.41   81.99\n');
        fprintf(fid,'Beam2RV= 278.79  -25.76   61.56  102.46   95.52\n');
        fprintf(fid,'Beam3RV= 273.64  -33.14   71.79  118.21   84.69\n');
        fprintf(fid,'Beam3LH= 286.6   -34.21   64.91  141.66   92.66\n');
        fprintf(fid,'Beam2LH= 264.78  -29.2    79.54  104.9    99.82\n');
        fprintf(fid,'Beam1LH= 266.27  -26.54   72.65  162.33   68.08\n');
        fprintf(fid,'Beam1RH= 261.47  -30.17   78.79   80.44   74.23\n');
        fprintf(fid,'Beam2RH= 268.69  -34.62   78.43   66.94  122.68\n');
        fprintf(fid,'Beam3RH= 266.93  -38.12   81.68   98.18  100.18\n\n');
    end
    fprintf(fid,'[adjustment]\n');
    fprintf(fid,['Date=',date,'\n']);
    fprintf(fid,'Description=Christoph Rudiger\n');
    fprintf(fid,['Comment=calibration using ',geb,' data (20',file_date(1:2),'.', ... 
        file_date(3:4),'.',file_date(5:6),')\n']);
    fprintf(fid,'Parameters=Offset,Slope\n');
    fprintf(fid,'Beam3LV= %15.12f %15.12f\n',off_vpol(1),slp_vpol(1));
    fprintf(fid,'Beam2LV= %15.12f %15.12f\n',off_vpol(2),slp_vpol(2));
    fprintf(fid,'Beam1LV= %15.12f %15.12f\n',off_vpol(3),slp_vpol(3));
    fprintf(fid,'Beam1RV= %15.12f %15.12f\n',off_vpol(4),slp_vpol(4));
    fprintf(fid,'Beam2RV= %15.12f %15.12f\n',off_vpol(5),slp_vpol(5));
    fprintf(fid,'Beam3RV= %15.12f %15.12f\n',off_vpol(6),slp_vpol(6));
    fprintf(fid,'Beam3LH= %15.12f %15.12f\n',off_hpol(1),slp_hpol(1));
    fprintf(fid,'Beam2LH= %15.12f %15.12f\n',off_hpol(2),slp_hpol(2));
    fprintf(fid,'Beam1LH= %15.12f %15.12f\n',off_hpol(3),slp_hpol(3));
    fprintf(fid,'Beam1RH= %15.12f %15.12f\n',off_hpol(4),slp_hpol(4));
    fprintf(fid,'Beam2RH= %15.12f %15.12f\n',off_hpol(5),slp_hpol(5));
    fprintf(fid,'Beam3RH= %15.12f %15.12f\n',off_hpol(6),slp_hpol(6));
    fclose(fid);
end
fclose('all')

cdpath=['cd ',fullfile(fpath,'Data',experi,'field',file_date,'PLMRcal')];
eval(cdpath)
copyfile([file_date,'_plmr.cal'],fullfile(fpath,'Data',experi,'field',file_date,'CUBE','plmr.cal'),'f');

end
