clear all
close all

%cd('/home/bandma01/Documents/WERAs/')
cd('/home/bandma01/hfr-data-processing/calibrations/')

site = 'PAO';
% site = 'PAB';

cd(site)

% system
for sys = 1:2

    if sys == 1
        systm = 'Helzel';
    elseif sys == 2
        systm = 'Flament';
    end

    clear J I E Phase calib_file calib_hr calib_jd calib_min calib_time calib_value
    clear calib_year d data delta_phase_obs delta_phase_th gps_data index_timecalib
    clear theta_pairs

    if strcmp(site,'PAO') == 1

        load('HelzelPAO_Track_29-OCT-14 183138.mat')

    elseif strcmp(site,'PAB') == 1

        if strcmp(systm,'Helzel') == 1
            load('HelzelPAB_Track_30-OCT-14 175446.mat')
        elseif strcmp(systm,'Flament') == 1
            load('FlamentPAB_Track_30-OCT-14 182008.mat')
        end

    end


    % phases and time of calibration
    if strcmp(site,'PAO') == 1
        cd('2014302')
    elseif strcmp(site,'PAB') == 1
        cd('2014303')
    end

    % list files
    clear d
    d = dir('*.mat');       % returns a structure with name, date, bytes, isdir, datenum
    Nfiles = length(d);     % nombre de fichiers a analyser


    % look at the names of files and extract information
    for c = 1:Nfiles
        calib_year(c) = str2double(d(c).name(1:4));
        calib_jd(c) = str2double(d(c).name(5:7));
        calib_hr(c) = str2double(d(c).name(8:9));
        calib_min(c) = str2double(d(c).name(10:11));

        load(d(c).name,'phase')

        Phase(:,c) = phase;
    end
    clear phase

    % for each time stamp (column) compute phase difference between 2
    % consecutive antennas (row)
    delta_phase_obs = modsym(diff(Phase,1),180);

    To = datenum(2014,01,01,0,0,0);

    % add 20 sec to center calibration time
    calib_time = To + calib_jd -1 + calib_hr./24 + calib_min/24/60 + 20/24/3600;
    % datevec(calib_time)

        clear calib_year calib_jd calib_hr calib_min

    if strcmp(site,'PAO') == 1
        % calib time 2min34sec late wrt GPS
        calib_time = calib_time + 2/24/60 + 34/24/3600;
    end
    disp(['calibrations started at ' datestr(calib_time(1)) ' and ended at ' datestr(calib_time(end))])

    % separate between calibrations made with Helzel's system and Flament's system
    % date of new calibration files
    % Helzel 2014-10-30 16:52 to 17:54
    % Flament 2014-10-30 17:55 to 18:20
    T1 = datenum(2014,10,30,17,55,00);

        clear J
        if strcmp(systm,'Helzel') == 1
            J = find(calib_time > T1);
        elseif strcmp(systm,'Flament') == 1
            J = find(calib_time < T1);
        end

        % remove unwanted data to avoid confusion!
%         calib_year(J)=[];
%         calib_jd(J)=[];
%         calib_hr(J)=[];
%         calib_min(J)=[];
        Phase(:,J)=[];
        delta_phase_obs(:,J)=[];
        calib_time(J)=[];
        d(J)=[];
        clear c Nfiles

    %+++++++++++++++++++++++++++++++++++++++++
    % for each calibration file  (series# 1)
    %+++++++++++++++++++++++++++++++++++++++++
    for calib_file = 1:length(calib_time)

        time_calib = calib_time(calib_file);

        % find corresponding file in gps data
        clear I
        I = find(data(:,1) > time_calib);
        % take closest point for now
        I = I(1);

        % datevec(data(I,1))
        % datevec(time_calib)

        index_timecalib(calib_file) = I;
    end

    % check
    % datevec(data(index_timecalib,1))

    % isolate gps calibration data
    gps_data = data(index_timecalib,:); % time lat lon

    figure(1)
    hold on
    h1=plot(gps_data(:,2),gps_data(:,3),'.');
    grid on
        if strcmp(systm,'Helzel') == 1
            set(h1,'marker','.','color','b')
        elseif strcmp(systm,'Flament') == 1
            set(h1,'marker','x','color','k')
        end

    % load antenna coordinates
    cd ..
    load([site '_Rx_coords.mat'])

    hold on
    plot(lon,lat,'x')
    grid on
    legend('gps','antennas','location','southeast')
    xlabel('longitude')
    ylabel('latitude')
    title(['site = ' site ])



    % compute theoretical phase values
    % wavelength
    lambda = 3E8/16.15E6;


    % initialize variables
    theta_pairs = nan(11,length(gps_data));

    % for each gps point
    for pt = 1:length(gps_data)

        % compute distance between GPS and antennas
        for ant = 1:12
            clear E W
            [E,W] = lonlat2km(gps_data(pt,2),gps_data(pt,3),lon(ant),lat(ant));
            ddist = sqrt(E.^2+W.^2);
            r(ant) = ddist*1000;    %meters
        end


        % for each pair of antennas
        for p = 1:11

            %delta_r
            delta_r(p) = r(p+1) - r(p);

            % theoretical phase
            delta_phase_th(p,pt) = modsym((delta_r(p)*360./lambda),180);


            % position of middle of pairs
            clear lon_mid lat_mid
            lon_mid = (lon(p+1)+lon(p))/2;
            lat_mid = (lat(p+1)+lat(p))/2;

            % distance between middle of pair and one of the antenna in the pair
            %         clear E W
            %         [E,W] = lonlat2km(lon_mid,lat_mid,lon(p),lat(p));
            %         ddist = sqrt(E.^2+W.^2);
            %         a(p) = ddist*1000;    %meters

            % distance between middle of pair and GPS pt
            clear ddsit phasAng
            [ddist,phasAng] = sw_dist([gps_data(pt,3) lat_mid],[gps_data(pt,2) lon_mid],'km');
            r_mid(p) = ddist*1000;    %meters
            zeta(p) = phasAng;


            % distance between pair of antennas
            clear ddsit phasAng
            [ddist,phasAng] = sw_dist([lat(p) lat(p+1)],[lon(p) lon(p+1)],'km');
            phi(p) = phasAng;


            % azimuth
            theta_pairs(p,pt) = modsym(zeta(p)-phi(p)+90,180);
            %rad2deg(theta_pairs)
        end
    end



    % for each pair of antenna
    for p = 1:11
        figure(p+1)
        h2=plot(theta_pairs(p,:),delta_phase_th(p,:),'x');
        hold on
        h3=plot(theta_pairs(p,:),-delta_phase_obs(p,:),'or');
        grid on
        xlabel('incidence (degrees)')
        ylabel('phase difference (degrees)')
        legend('theoretical','observed','location','best')
        if strcmp(systm,'Helzel') == 1
            set(h2,'color','b')
            set(h3,'color','r','marker','o')
        elseif strcmp(systm,'Flament') == 1
            set(h2,'color','k','marker','*')
            set(h3,'color','r','marker','sq','markerfacecolor','r')
        end
        title([site ' - antennas ' num2str(p) '-' num2str(p+1)])


        % compute calibration value for antennas pairs
        calib_value(p) = nanmean(modsym(-delta_phase_obs(p,:)-delta_phase_th(p,:),180),2);


         % compute calibration value for each antenna
        phi_ant(p+1,sys) = sum(calib_value);

    end

    P = ones(16,3);
    P(13:16,:) = 0;
    P(1:12,3) = phi_ant(:,1);


    % save file (we only keep results from Helzel
    if strcmp(site,'PAO') == 1
        dlmwrite('calibration_pao_29OCT2014.wera',P,'delimiter','\t','precision','%.7f')
    elseif strcmp(site,'PAB') == 1
        dlmwrite('calibration_pab_30OCT2014.wera',P,'delimiter','\t','precision','%.7f')
    end


    figure(13)
    hold on
    h4 = plot(calib_value,'.-');
    grid on
    xlabel('pairs of antennas')
    ylabel(['phase differences between theoretical',10 ,'and observed values (degrees)'])
        if strcmp(systm,'Helzel') == 1
            set(h4,'color','b')
            legend('Helzel')
        elseif strcmp(systm,'Flament') == 1
            set(h4,'color','k')
            legend('Helzel','Flament','location','southeast')
        end
    title(site)

    figure(14)
    hold on
    h5 = plot([1:12],phi_ant(:,sys),'r');
    grid on
    if strcmp(systm,'Helzel') == 1
        set(h5,'marker','o')
        legend('Helzel')
    elseif strcmp(systm,'Flament') == 1
        set(h5,'marker','sq','markerfacecolor','r')
        legend('Helzel','Flament','location','southeast')
    end
    xlabel('antenna number')
    ylabel('calibration phase to be applied (degrees)')
    title(site)

    if strcmp(site,'PAO') == 1
       break
    end

end

% add current calibration phases
if strcmp(site,'PAO') == 1
    load calibration_pao.wera
    figure(14)
    plot([1:12],calibration_pao(1:12,3),'b.-')
    legend('Helzel','original calibration','location','northwest')
elseif strcmp(site,'PAB') == 1
    load calibration_pab.wera
    figure(14)
    plot([1:12],calibration_pab(1:12,3),'b.-')
    legend('Helzel','Flament','original calibration','location','southeast')
end






% save figures
save_flag = 1;

figdir = ('/home/bandma01/Dropbox/UQAR/MEOPAR_Marion/figures/WERAcalibration');
cd(figdir)

fig_hdl = figure(1);
if save_flag
    fname = ['configuration_' site '_calibrations'];
    fig_w = 8;%cm
    fig_h = 5;%cm
    set(fig_hdl,'PaperUnits','centimeters')
    set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])
    set(fig_hdl, 'PaperSize', [fig_w fig_h]);
    print(fig_hdl,'-dpng','-r500',fname)
    %crop([fname '.png'])
end




for i = 2:12
    fig_hdl = figure(i);
    if save_flag
        fname = ['Phases_' site '_pair' num2str(i-1) '-' num2str(i)];
        fig_w = 8;%cm
        fig_h = 5;%cm
        set(fig_hdl,'PaperUnits','centimeters')
        set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])
        set(fig_hdl, 'PaperSize', [fig_w fig_h]);
        print(fig_hdl,'-dpng','-r500',fname)
        %         crop([fname '.png'])
    end
end



fig_hdl = figure(13);
if save_flag
    fname = ['PhaseDifference_' site '_calibrations'];
    fig_w = 8;%cm
    fig_h = 5;%cm
    set(fig_hdl,'PaperUnits','centimeters')
    set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])
    set(fig_hdl, 'PaperSize', [fig_w fig_h]);
    print(fig_hdl,'-dpng','-r500',fname)
    %         crop([fname '.png'])
end

fig_hdl = figure(14);
if save_flag
    fname = ['PhaseCalibration_' site];
    fig_w = 8;%cm
    fig_h = 5;%cm
    set(fig_hdl,'PaperUnits','centimeters')
    set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])
    set(fig_hdl, 'PaperSize', [fig_w fig_h]);
    print(fig_hdl,'-dpng','-r500',fname)
    %         crop([fname '.png'])
end

calibrations started at 29-Oct-2014 17:48:54 and ended at 29-Oct-2014 18:25:54

Published with MATLAB® 7.14