Commit e9cb9445f568915ba954ee3d7397d43080eef063

Authored by Marion Bandet
1 parent 7e281bfb
Exists in master

update calibrations PAB 2016

Showing 99 changed files with 1340 additions and 252 deletions   Show diff stats
build/doctrees/calibPAB2014.doctree
No preview for this file type
build/doctrees/calibPAB2016.doctree 0 → 100644
No preview for this file type
build/doctrees/calibPAO2014.doctree
No preview for this file type
build/doctrees/calibration.doctree
No preview for this file type
build/doctrees/environment.pickle
No preview for this file type
build/doctrees/index.doctree
No preview for this file type
build/html/RenameFiles.sh
... ... @@ -14,8 +14,8 @@ set -x
14 14 # mv "$f" "$(echo "$f" | sed -e 's/_/_STLE_/' - )"
15 15 #done
16 16  
17   -for f in WERAcalibration_*.png
  17 +for f in WERAcalibration_2016*.png
18 18 do
19 19 ### mv "$f" "$(echo "$f" | sed -e 's/[:]//' - )"
20   - mv "$f" "$(echo "$f" | sed -e 's/WERAcalibration_2014/WERAcalibrationPAO_2014/' - )"
  20 + mv "$f" "$(echo "$f" | sed -e 's/WERAcalibration_2016/WERAcalibrationPAB_2016/' - )"
21 21 done
... ...
build/html/WERAcalibrationPAB_2016.html 0 → 100644
... ... @@ -0,0 +1,1001 @@
  1 +
  2 +<!DOCTYPE html
  3 + PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
  4 +<html><head>
  5 + <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
  6 + <!--
  7 +This HTML was auto-generated from MATLAB code.
  8 +To make changes, update the MATLAB code and republish this document.
  9 + --><title>WERAcalibration_2016</title><meta name="generator" content="MATLAB 7.14"><link rel="schema.DC" href="http://purl.org/dc/elements/1.1/"><meta name="DC.date" content="2016-07-22"><meta name="DC.source" content="WERAcalibration_2016.m"><style type="text/css">
  10 +html,body,div,span,applet,object,iframe,h1,h2,h3,h4,h5,h6,p,blockquote,pre,a,abbr,acronym,address,big,cite,code,del,dfn,em,font,img,ins,kbd,q,s,samp,small,strike,strong,sub,sup,tt,var,b,u,i,center,dl,dt,dd,ol,ul,li,fieldset,form,label,legend,table,caption,tbody,tfoot,thead,tr,th,td{margin:0;padding:0;border:0;outline:0;font-size:100%;vertical-align:baseline;background:transparent}body{line-height:1}ol,ul{list-style:none}blockquote,q{quotes:none}blockquote:before,blockquote:after,q:before,q:after{content:'';content:none}:focus{outine:0}ins{text-decoration:none}del{text-decoration:line-through}table{border-collapse:collapse;border-spacing:0}
  11 +
  12 +html { min-height:100%; margin-bottom:1px; }
  13 +html body { height:100%; margin:0px; font-family:Arial, Helvetica, sans-serif; font-size:10px; color:#000; line-height:140%; background:#fff none; overflow-y:scroll; }
  14 +html body td { vertical-align:top; text-align:left; }
  15 +
  16 +h1 { padding:0px; margin:0px 0px 25px; font-family:Arial, Helvetica, sans-serif; font-size:1.5em; color:#d55000; line-height:100%; font-weight:normal; }
  17 +h2 { padding:0px; margin:0px 0px 8px; font-family:Arial, Helvetica, sans-serif; font-size:1.2em; color:#000; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }
  18 +h3 { padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:1.1em; color:#000; font-weight:bold; line-height:140%; }
  19 +
  20 +a { color:#005fce; text-decoration:none; }
  21 +a:hover { color:#005fce; text-decoration:underline; }
  22 +a:visited { color:#004aa0; text-decoration:none; }
  23 +
  24 +p { padding:0px; margin:0px 0px 20px; }
  25 +img { padding:0px; margin:0px 0px 20px; border:none; }
  26 +p img, pre img, tt img, li img { margin-bottom:0px; }
  27 +
  28 +ul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }
  29 +ul li { padding:0px; margin:0px 0px 7px 0px; }
  30 +ul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }
  31 +ul li ol li { list-style:decimal; }
  32 +ol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }
  33 +ol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }
  34 +ol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }
  35 +ol li ol li { list-style-type:lower-alpha; }
  36 +ol li ul { padding-top:7px; }
  37 +ol li ul li { list-style:square; }
  38 +
  39 +.content { font-size:1.2em; line-height:140%; padding: 20px; }
  40 +
  41 +pre, tt, code { font-size:12px; }
  42 +pre { margin:0px 0px 20px; }
  43 +pre.error { color:red; }
  44 +pre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }
  45 +pre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }
  46 +
  47 +@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }
  48 +
  49 +span.keyword { color:#0000FF }
  50 +span.comment { color:#228B22 }
  51 +span.string { color:#A020F0 }
  52 +span.untermstring { color:#B20000 }
  53 +span.syscmd { color:#B28C00 }
  54 +
  55 +.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }
  56 +.footer p { margin:0px; }
  57 +
  58 + </style></head><body><div class="content"><h2>Contents</h2><div><ul><li><a href="#3">separate between calibrations made walking and with the boat</a></li><li><a href="#4">compute theoretical phase values</a></li></ul></div><pre class="codeinput">clear <span class="string">all</span>
  59 +close <span class="string">all</span>
  60 +
  61 +<span class="comment">% calibrations made with the Helzel system</span>
  62 +
  63 +
  64 +<span class="comment">%cd('/home/bandma01/Documents/WERAs/')</span>
  65 +cd(<span class="string">'/home/bandma01/hfr-data-processing/calibrations/'</span>)
  66 +<span class="comment">% datadir = '/run/media/bandma01/SeagateWireless/MEOPAR/DATA/HFRs/WERAs/';</span>
  67 +<span class="comment">% cd(datadir)</span>
  68 +
  69 +<span class="comment">% site = 'PAO';</span>
  70 +site = <span class="string">'PAB'</span>;
  71 +
  72 +cd(site)
  73 +
  74 +<span class="comment">% calibrations: WALK or BOAT1 or BOAT2</span>
  75 +<span class="keyword">for</span> sys = 1:3
  76 +</pre><pre class="codeinput"> <span class="keyword">if</span> sys == 1
  77 + systm = <span class="string">'WALK'</span>;
  78 + <span class="keyword">elseif</span> sys == 2
  79 + systm = <span class="string">'BOAT1'</span>;
  80 + <span class="keyword">elseif</span> sys == 3
  81 + systm = <span class="string">'BOAT2'</span>;
  82 + <span class="keyword">end</span>
  83 +
  84 + clear <span class="string">J</span> <span class="string">I</span> <span class="string">E</span> <span class="string">Phase</span> <span class="string">calib_file</span> <span class="string">calib_hr</span> <span class="string">calib_jd</span> <span class="string">calib_min</span> <span class="string">calib_time</span> <span class="string">calib_value</span>
  85 + clear <span class="string">calib_year</span> <span class="string">d</span> <span class="string">data</span> <span class="string">delta_phase_obs</span> <span class="string">delta_phase_th</span> <span class="string">gps_data</span> <span class="string">index_timecalib</span>
  86 + clear <span class="string">theta_pairs</span>
  87 +
  88 + <span class="keyword">if</span> strcmp(site,<span class="string">'PAO'</span>) == 1
  89 +
  90 + <span class="comment">% load('HelzelPAO_Track_29-OCT-14 183138.mat')</span>
  91 +
  92 + <span class="keyword">elseif</span> strcmp(site,<span class="string">'PAB'</span>) == 1
  93 +
  94 + <span class="keyword">if</span> strcmp(systm,<span class="string">'WALK'</span>) == 1
  95 + load(<span class="string">'Track_05-JUL-16 150101.mat'</span>)
  96 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT1'</span>) == 1
  97 + load(<span class="string">'FirstCircle_BOAT_Track_05-JUL-16.mat'</span>)
  98 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT2'</span>) == 1
  99 + load(<span class="string">'SecondCircle_BOAT_Track_05-JUL-16.mat'</span>)
  100 + <span class="keyword">end</span>
  101 +
  102 + <span class="keyword">end</span>
  103 +
  104 +
  105 + <span class="comment">% phases and time of calibration</span>
  106 + <span class="keyword">if</span> strcmp(site,<span class="string">'PAO'</span>) == 1
  107 + <span class="comment">%cd('2014302')</span>
  108 + <span class="keyword">elseif</span> strcmp(site,<span class="string">'PAB'</span>) == 1
  109 + cd(<span class="string">'2016187'</span>)
  110 + <span class="keyword">end</span>
  111 +
  112 + <span class="comment">% list files</span>
  113 + clear <span class="string">d</span>
  114 + d = dir(<span class="string">'*.mat'</span>); <span class="comment">% returns a structure with name, date, bytes, isdir, datenum</span>
  115 + Nfiles = length(d); <span class="comment">% nombre de fichiers a analyser</span>
  116 +
  117 +
  118 + <span class="comment">% look at the names of files and extract information</span>
  119 + <span class="keyword">for</span> c = 1:Nfiles
  120 + calib_year(c) = str2double(d(c).name(1:4));
  121 + calib_jd(c) = str2double(d(c).name(5:7));
  122 + calib_hr(c) = str2double(d(c).name(8:9));
  123 + calib_min(c) = str2double(d(c).name(10:11));
  124 +
  125 + load(d(c).name,<span class="string">'phase'</span>)
  126 +
  127 + Phase(:,c) = phase;
  128 + <span class="keyword">end</span>
  129 + clear <span class="string">phase</span>
  130 +
  131 + <span class="comment">% for each time stamp (column) compute phase difference between 2</span>
  132 + <span class="comment">% consecutive antennas (row)</span>
  133 + delta_phase_obs = modsym(diff(Phase,1),180);
  134 +
  135 + To = datenum(2016,01,01,0,0,0);
  136 +
  137 +<span class="comment">% % add 20 sec to center calibration time</span>
  138 +<span class="comment">% calib_time = To + calib_jd -1 + calib_hr./24 + calib_min/24/60 + 20/24/3600;</span>
  139 +<span class="comment">% % datevec(calib_time)</span>
  140 +<span class="comment">%</span>
  141 +<span class="comment">% clear calib_year calib_jd calib_hr calib_min</span>
  142 +
  143 + calib_time = To + calib_jd -1 + calib_hr./24 + calib_min/24/60;
  144 +
  145 + <span class="keyword">if</span> strcmp(site,<span class="string">'PAO'</span>) == 1
  146 + <span class="comment">% calib time 2min34sec late wrt GPS</span>
  147 + calib_time = calib_time + 2/24/60 + 34/24/3600;
  148 + <span class="keyword">end</span>
  149 + disp([<span class="string">'the acquisition of calibration started at '</span> datestr(calib_time(1)) <span class="string">' and ended up at '</span> datestr(calib_time(end))])
  150 +</pre><pre class="codeoutput">the acquisition of calibration started at 05-Jul-2016 00:13:00 and ended up at 05-Jul-2016 20:11:00
  151 +</pre><pre class="codeoutput">the acquisition of calibration started at 05-Jul-2016 00:13:00 and ended up at 05-Jul-2016 20:11:00
  152 +</pre><pre class="codeoutput">the acquisition of calibration started at 05-Jul-2016 00:13:00 and ended up at 05-Jul-2016 20:11:00
  153 +</pre><h2>separate between calibrations made walking and with the boat<a name="3"></a></h2><p>date of new calibration files walk 2016-07-05 13:53 and 14:58 UTC boat 1st circle 2016-07-05 18:20 and 19:18 UTC boat 2nd circle 2016-07-05 19:20 and 20:11 UTC</p><pre class="codeinput"> <span class="keyword">if</span> strcmp(systm,<span class="string">'WALK'</span>) == 1
  154 + T1 = datenum(2016,07,05,13,53,00);
  155 + T2 = datenum(2016,07,05,14,58,00);
  156 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT1'</span>) == 1
  157 + T1 = datenum(2016,07,05,18,20,00);
  158 + T2 = datenum(2016,07,05,19,18,00);
  159 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT2'</span>) == 1
  160 + T1 = datenum(2016,07,05,19,20,00);
  161 + T2 = datenum(2016,07,05,20,11,00);
  162 + <span class="keyword">end</span>
  163 +
  164 + clear <span class="string">J</span> <span class="string">K</span>
  165 + J = find(calib_time &lt; T1);
  166 + K = find(calib_time &gt; T2);
  167 + J = [J K]; clear <span class="string">K</span>
  168 +
  169 + <span class="comment">% remove unwanted data to avoid confusion!</span>
  170 +<span class="comment">% calib_year(J)=[];</span>
  171 +<span class="comment">% calib_jd(J)=[];</span>
  172 +<span class="comment">% calib_hr(J)=[];</span>
  173 +<span class="comment">% calib_min(J)=[];</span>
  174 + Phase(:,J)=[];
  175 + delta_phase_obs(:,J)=[];
  176 + calib_time(J)=[];
  177 + d(J)=[];
  178 + clear <span class="string">c</span> <span class="string">Nfiles</span>
  179 +
  180 + <span class="comment">%+++++++++++++++++++++++++++++++++++++++++</span>
  181 + <span class="comment">% for each calibration file (series# 1)</span>
  182 + <span class="comment">%+++++++++++++++++++++++++++++++++++++++++</span>
  183 + <span class="keyword">for</span> calib_file = 1:length(calib_time)
  184 +
  185 + time_calib = calib_time(calib_file);
  186 + <span class="comment">% datevec(time_calib)</span>
  187 +
  188 + <span class="comment">% find corresponding file in gps data</span>
  189 + clear <span class="string">I</span>
  190 + I = find(time_calib &lt;= data(:,1) );
  191 + <span class="comment">% take closest point for now</span>
  192 + I = I(1);
  193 +
  194 + <span class="comment">% datevec(data(I,1))</span>
  195 + <span class="comment">% datevec(time_calib)</span>
  196 +
  197 + index_timecalib(calib_file) = I;
  198 + <span class="keyword">end</span>
  199 +
  200 + <span class="comment">% check</span>
  201 + <span class="comment">%datevec(data(index_timecalib,1))</span>
  202 +
  203 + <span class="comment">% isolate gps calibration data</span>
  204 + gps_data = data(index_timecalib,:); <span class="comment">% time lat lon</span>
  205 +
  206 + figure(1)
  207 + hold <span class="string">on</span>
  208 + <span class="keyword">if</span> sys == 1
  209 + h1(sys) = plot(gps_data(:,2),gps_data(:,3),<span class="string">'.-'</span>,<span class="string">'displayName'</span>,<span class="string">'gps WALK'</span>);
  210 + <span class="keyword">elseif</span> sys == 2
  211 + h1(sys) = plot(gps_data(:,2),gps_data(:,3),<span class="string">'.-'</span>,<span class="string">'displayName'</span>,<span class="string">'gps BOAT circle 1'</span>);
  212 + <span class="keyword">elseif</span> sys == 3
  213 + h1(sys) = plot(gps_data(:,2),gps_data(:,3),<span class="string">'.-'</span>,<span class="string">'displayName'</span>,<span class="string">'gps BOAT circle 2'</span>);
  214 + <span class="keyword">end</span>
  215 +
  216 + grid <span class="string">on</span>
  217 +
  218 + <span class="keyword">if</span> strcmp(systm,<span class="string">'WALK'</span>) == 1
  219 + set(h1(1),<span class="string">'marker'</span>,<span class="string">'x'</span>,<span class="string">'color'</span>,<span class="string">'b'</span>)
  220 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT1'</span>) == 1
  221 + set(h1(2),<span class="string">'marker'</span>,<span class="string">'.'</span>,<span class="string">'color'</span>,[0 0.5 0])
  222 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT2'</span>) == 1
  223 + set(h1(3),<span class="string">'marker'</span>,<span class="string">'*'</span>,<span class="string">'color'</span>,[0.6 0.2 0])
  224 + <span class="keyword">end</span>
  225 +
  226 + <span class="comment">% load antenna coordinates</span>
  227 + <span class="comment">%</span>
  228 + <span class="comment">%load([site '_Rx_coords.mat'])</span>
  229 +
  230 +
  231 + <span class="keyword">if</span> strcmp(site,<span class="string">'PAO'</span>) == 1
  232 + cd(<span class="string">'/home/bandma01/hfr-data-processing/calibrations/'</span>)
  233 + load([<span class="string">'Survey_2014_'</span> site])
  234 + lon = PAO.lon;
  235 + lat = PAO.lat;
  236 + <span class="keyword">elseif</span> strcmp(site,<span class="string">'PAB'</span>) == 1
  237 + cd <span class="string">..</span>
  238 + load([site <span class="string">'_Rx_coords.mat'</span>])
  239 + <span class="keyword">end</span>
  240 +
  241 + hold <span class="string">on</span>
  242 + h2b(sys) = plot(lon,lat,<span class="string">'dk'</span>,<span class="string">'displayName'</span>,<span class="string">'antennas positions'</span>);
  243 + grid <span class="string">on</span>
  244 + <span class="keyword">if</span> strcmp(systm,<span class="string">'WALK'</span>) == 1
  245 + <span class="comment">%legend('gps Walk','antennas positions','location','best')</span>
  246 + legend([h2b h1],<span class="string">'location'</span>,<span class="string">'best'</span>)
  247 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT1'</span>) == 1
  248 + legend([h2b(1) h1],<span class="string">'location'</span>,<span class="string">'best'</span>)
  249 + <span class="comment">% legend('gps Walk','antennas positions','gps Boat circle 1','antennas positions','location','best')</span>
  250 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT2'</span>) == 1
  251 + legend([h2b(1) h1],<span class="string">'location'</span>,<span class="string">'best'</span>)
  252 + <span class="comment">% legend('gps Walk','antennas positions','gps Boat circle 1','antennas positions','gps boat circle 2','antennas positions','location','best')</span>
  253 + <span class="keyword">end</span>
  254 + xlabel(<span class="string">'longitude'</span>)
  255 + ylabel(<span class="string">'latitude'</span>)
  256 + title([<span class="string">'site = '</span> site ])
  257 +
  258 +<span class="comment">% continue</span>
  259 +</pre><img vspace="5" hspace="5" src="WERAcalibration_2016_01.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_15.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_29.png" alt=""> <h2>compute theoretical phase values<a name="4"></a></h2><p>wavelength</p><pre class="codeinput"> lambda = 3E8/16.15E6;
  260 +
  261 +
  262 + <span class="comment">% initialize variables</span>
  263 + theta_pairs = nan(11,length(gps_data));
  264 +
  265 + <span class="comment">% for each gps point</span>
  266 + <span class="keyword">for</span> pt = 1:length(gps_data)
  267 +
  268 + <span class="comment">% compute distance between GPS and antennas</span>
  269 + <span class="keyword">for</span> ant = 1:12
  270 + clear <span class="string">E</span> <span class="string">W</span>
  271 + [E,W] = lonlat2km(gps_data(pt,2),gps_data(pt,3),lon(ant),lat(ant));
  272 + ddist = sqrt(E.^2+W.^2);
  273 + r(ant) = ddist*1000; <span class="comment">%meters</span>
  274 + <span class="keyword">end</span>
  275 +
  276 +
  277 + <span class="comment">% for each pair of antennas</span>
  278 + <span class="keyword">for</span> p = 1:11
  279 +
  280 + <span class="comment">%delta_r</span>
  281 + delta_r(p) = r(p+1) - r(p);
  282 +
  283 + <span class="comment">% theoretical phase</span>
  284 + delta_phase_th(p,pt) = modsym((delta_r(p)*360./lambda),180);
  285 +
  286 +
  287 + <span class="comment">% position of middle of pairs</span>
  288 + clear <span class="string">lon_mid</span> <span class="string">lat_mid</span>
  289 + lon_mid = (lon(p+1)+lon(p))/2;
  290 + lat_mid = (lat(p+1)+lat(p))/2;
  291 +
  292 + <span class="comment">% distance between middle of pair and one of the antenna in the pair</span>
  293 + <span class="comment">% clear E W</span>
  294 + <span class="comment">% [E,W] = lonlat2km(lon_mid,lat_mid,lon(p),lat(p));</span>
  295 + <span class="comment">% ddist = sqrt(E.^2+W.^2);</span>
  296 + <span class="comment">% a(p) = ddist*1000; %meters</span>
  297 +
  298 + <span class="comment">% distance between middle of pair and GPS pt</span>
  299 + clear <span class="string">ddsit</span> <span class="string">phasAng</span>
  300 + [ddist,phasAng] = sw_dist([gps_data(pt,3) lat_mid],[gps_data(pt,2) lon_mid],<span class="string">'km'</span>);
  301 + r_mid(p) = ddist*1000; <span class="comment">%meters</span>
  302 + zeta(p) = phasAng;
  303 +
  304 +
  305 + <span class="comment">% distance between pair of antennas</span>
  306 + clear <span class="string">ddsit</span> <span class="string">phasAng</span>
  307 + [ddist,phasAng] = sw_dist([lat(p) lat(p+1)],[lon(p) lon(p+1)],<span class="string">'km'</span>);
  308 + phi(p) = phasAng;
  309 +
  310 +
  311 + <span class="comment">% azimuth</span>
  312 + theta_pairs(p,pt) = modsym(zeta(p)-phi(p)+90,180);
  313 + <span class="comment">%rad2deg(theta_pairs)</span>
  314 + <span class="keyword">end</span>
  315 + <span class="keyword">end</span>
  316 +
  317 +
  318 +
  319 + <span class="comment">% for each pair of antenna</span>
  320 + <span class="keyword">for</span> p = 1:11
  321 +
  322 + figure(p+1)
  323 +
  324 + h2(sys,p) = plot(theta_pairs(p,:),delta_phase_th(p,:),<span class="string">'-'</span>,<span class="string">'displayName'</span>,<span class="string">'theoretical'</span>);
  325 + hold <span class="string">on</span>
  326 + <span class="keyword">if</span> sys == 1
  327 + h3(sys,p) = plot(theta_pairs(p,:),-delta_phase_obs(p,:),<span class="string">'o'</span>,<span class="string">'displayName'</span>,<span class="string">'WALK calib'</span>);
  328 + <span class="keyword">elseif</span> sys == 2
  329 + h3(sys,p) = plot(theta_pairs(p,:),-delta_phase_obs(p,:),<span class="string">'o'</span>,<span class="string">'displayName'</span>,<span class="string">'BOAT circle 1 calib'</span>);
  330 + <span class="keyword">elseif</span> sys == 3
  331 + h3(sys,p) = plot(theta_pairs(p,:),-delta_phase_obs(p,:),<span class="string">'o'</span>,<span class="string">'displayName'</span>,<span class="string">'BOAT circle 2 calib'</span>);
  332 + <span class="keyword">end</span>
  333 + grid <span class="string">on</span>
  334 + xlabel(<span class="string">'incidence (degrees)'</span>)
  335 + ylabel(<span class="string">'phase difference (degrees)'</span>)
  336 +
  337 + <span class="keyword">if</span> strcmp(systm,<span class="string">'WALK'</span>) == 1
  338 + set(h2(sys,p),<span class="string">'color'</span>,<span class="string">'b'</span>)
  339 + set(h3(sys,p),<span class="string">'color'</span>,<span class="string">'r'</span>,<span class="string">'marker'</span>,<span class="string">'o'</span>)
  340 + legend([h2(1,p) h3(1,p)],<span class="string">'location'</span>,<span class="string">'best'</span>)
  341 + <span class="comment">%legend('theoretical','Walk','location','best')</span>
  342 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT1'</span>) == 1
  343 + set(h2(sys,p),<span class="string">'color'</span>,<span class="string">'b'</span>)
  344 + set(h3(sys,p),<span class="string">'color'</span>,<span class="string">'r'</span>,<span class="string">'marker'</span>,<span class="string">'sq'</span>,<span class="string">'markerfacecolor'</span>,<span class="string">'r'</span>)
  345 + <span class="comment">%legend('theoretical','Boat circle 1','location','best')</span>
  346 + legend([h2(1,p) h3(1,p) h3(2,p)],<span class="string">'location'</span>,<span class="string">'best'</span>)
  347 + <span class="keyword">elseif</span> strcmp(systm,<span class="string">'BOAT2'</span>) == 1
  348 + set(h2(sys,p),<span class="string">'color'</span>,<span class="string">'b'</span>)
  349 + set(h3(sys,p),<span class="string">'color'</span>,[0 0.5 0],<span class="string">'marker'</span>,<span class="string">'sq'</span>,<span class="string">'markerfacecolor'</span>,[0 0.5 0])
  350 + <span class="comment">%legend('theoretical','Boat circle 2','location','best')</span>
  351 + legend([h2(1,p) h3(1,p) h3(2,p) h3(3,p)],<span class="string">'location'</span>,<span class="string">'best'</span>)
  352 + <span class="keyword">end</span>
  353 + title([site <span class="string">' - antennas '</span> num2str(p) <span class="string">'-'</span> num2str(p+1)])
  354 +
  355 +
  356 + <span class="comment">% compute calibration value for antennas pairs</span>
  357 + <span class="comment">% calib_value(p) = nanmean(modsym(-delta_phase_obs(p,:)-delta_phase_th(p,:),180),2);</span>
  358 + calib_value(p) = nanmean(modsym(delta_phase_th(p,:)+delta_phase_obs(p,:),180),2);
  359 +
  360 +
  361 + <span class="comment">% compute calibration value for each antenna</span>
  362 + phi_ant(p+1,sys) = sum(calib_value);
  363 +
  364 + <span class="keyword">end</span>
  365 +
  366 + <span class="comment">%continue</span>
  367 +
  368 + P = ones(16,3);
  369 + P(13:16,:) = 0;
  370 + P(1:12,3) = phi_ant(:,1);
  371 +
  372 +
  373 + <span class="comment">% save file (we only keep results from Helzel</span>
  374 +<span class="comment">% datasav = [datadir site];</span>
  375 +<span class="comment">% cd(datasav)</span>
  376 +<span class="comment">% if strcmp(site,'PAO') == 1</span>
  377 +<span class="comment">% dlmwrite('calibration_pao_29OCT2014.wera',P,'delimiter','\t','precision','%.7f')</span>
  378 +<span class="comment">% elseif strcmp(site,'PAB') == 1</span>
  379 +<span class="comment">% dlmwrite('calibration_pab_30OCT2014.wera',P,'delimiter','\t','precision','%.7f')</span>
  380 +<span class="comment">% end</span>
  381 +
  382 +
  383 + figure(13)
  384 + hold <span class="string">on</span>
  385 + <span class="keyword">if</span> sys == 1
  386 + h4(sys) = plot(calib_value,<span class="string">'.-'</span>,<span class="string">'displayname'</span>,<span class="string">'WALK'</span>,<span class="string">'color'</span>,<span class="string">'r'</span>,<span class="string">'marker'</span>,<span class="string">'o'</span>);
  387 + <span class="keyword">elseif</span> sys == 2
  388 + h4(sys) = plot(calib_value,<span class="string">'.-'</span>,<span class="string">'displayname'</span>,<span class="string">'BOAT circle 1'</span>,<span class="string">'color'</span>,<span class="string">'k'</span>,<span class="string">'marker'</span>,<span class="string">'sq'</span>);
  389 + <span class="keyword">elseif</span> sys == 3
  390 + h4(sys) = plot(calib_value,<span class="string">'.-'</span>,<span class="string">'displayname'</span>,<span class="string">'BOAT circle 2'</span>,<span class="string">'color'</span>,[0 0.5 0],<span class="string">'marker'</span>,<span class="string">'d'</span>);
  391 + <span class="keyword">end</span>
  392 + legend(h4)
  393 + grid <span class="string">on</span>
  394 + xlabel(<span class="string">'pairs of antennas'</span>)
  395 + ylabel([<span class="string">'averaged phase differences between theoretical'</span>,10 ,<span class="string">'and observed values (degrees)'</span>])
  396 + title(site)
  397 +
  398 +
  399 +
  400 + figure(14)
  401 + hold <span class="string">on</span>
  402 + <span class="keyword">if</span> sys == 1
  403 + h5(sys) = plot([1:12],phi_ant(:,sys),<span class="string">'displayname'</span>,<span class="string">'WALK'</span>,<span class="string">'color'</span>,<span class="string">'r'</span>,<span class="string">'marker'</span>,<span class="string">'o'</span>);
  404 + <span class="keyword">elseif</span> sys == 2
  405 + h5(sys) = plot([1:12],phi_ant(:,sys),<span class="string">'displayname'</span>,<span class="string">'BOAT circle 1'</span>,<span class="string">'color'</span>,<span class="string">'k'</span>,<span class="string">'marker'</span>,<span class="string">'sq'</span>);
  406 + <span class="keyword">elseif</span> sys == 3
  407 + h5(sys) = plot([1:12],phi_ant(:,sys),<span class="string">'displayname'</span>,<span class="string">'BOAT circle 2'</span>,<span class="string">'color'</span>,[0 0.5 0],<span class="string">'marker'</span>,<span class="string">'d'</span>);
  408 + <span class="keyword">end</span>
  409 + grid <span class="string">on</span>
  410 + legend(h5)
  411 + xlabel(<span class="string">'antenna number'</span>)
  412 + ylabel(<span class="string">'calibration phase to be applied (degrees)'</span>)
  413 + title(site)
  414 +
  415 + <span class="keyword">if</span> strcmp(site,<span class="string">'PAO'</span>) == 1
  416 + <span class="keyword">break</span>
  417 + <span class="keyword">end</span>
  418 +</pre><img vspace="5" hspace="5" src="WERAcalibration_2016_02.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_03.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_04.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_05.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_06.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_07.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_08.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_09.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_10.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_11.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_12.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_13.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_14.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_16.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_17.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_18.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_19.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_20.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_21.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_22.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_23.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_24.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_25.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_26.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_27.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_28.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_30.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_31.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_32.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_33.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_34.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_35.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_36.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_37.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_38.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_39.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_40.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_41.png" alt=""> <img vspace="5" hspace="5" src="WERAcalibration_2016_42.png" alt=""> <pre class="codeinput"><span class="keyword">end</span>
  419 +
  420 +
  421 +<span class="comment">% add current calibration phases</span>
  422 +<span class="keyword">if</span> strcmp(site,<span class="string">'PAO'</span>) == 1
  423 + load <span class="string">calibration_pao.wera</span>
  424 + figure(14)
  425 + hold <span class="string">on</span>
  426 + plot([1:12],calibration_pao(1:12,3),<span class="string">'b.-'</span>)
  427 + legend(<span class="string">'Helzel'</span>,<span class="string">'original calibration'</span>,<span class="string">'location'</span>,<span class="string">'northwest'</span>)
  428 +<span class="keyword">elseif</span> strcmp(site,<span class="string">'PAB'</span>) == 1
  429 + load <span class="string">calibration_pab.wera</span>
  430 + figure(14)
  431 + hold <span class="string">on</span>
  432 + h6 = plot([1:12],calibration_pab(1:12,3),<span class="string">'b.-'</span>,<span class="string">'displayname'</span>,<span class="string">'primary calibrations'</span>)
  433 +
  434 + load <span class="string">calibration_pab_30OCT2014.wera</span>
  435 + figure(14)
  436 + hold <span class="string">on</span>
  437 + h7 = plot([1:12],calibration_pab_30OCT2014(1:12,3),<span class="string">'m.-'</span>,<span class="string">'displayname'</span>,<span class="string">'calibrations 30OCT2014'</span>)
  438 + legend([h5 h6 h7],<span class="string">'location'</span>,<span class="string">'best'</span>)
  439 +
  440 +<span class="keyword">end</span>
  441 +
  442 +
  443 +
  444 +
  445 +
  446 +
  447 +<span class="comment">% save figures</span>
  448 +save_flag = 0;
  449 +
  450 +<span class="comment">%figdir = ('/home/bandma01/Dropbox/UQAR/MEOPAR_Marion/figures/WERAcalibration');</span>
  451 +figdir = ([<span class="string">'/home/bandma01/hfr-data-processing/calibrations/'</span> site]);
  452 +
  453 +cd(figdir)
  454 +
  455 +fig_hdl = figure(1);
  456 +
  457 +<span class="comment">% Adjust Font and Axes Properties</span>
  458 +set( gca , <span class="keyword">...</span>
  459 + <span class="string">'FontName'</span> , <span class="string">'Helvetica'</span>);
  460 +
  461 +
  462 +<span class="keyword">if</span> save_flag
  463 + fname = [<span class="string">'configuration_'</span> site <span class="string">'_calibrations'</span>];
  464 + <span class="comment">%fig_w = 8;%cm</span>
  465 + <span class="comment">%fig_h = 5;%cm</span>
  466 + <span class="comment">%set(fig_hdl,'PaperUnits','centimeters')</span>
  467 + <span class="comment">%set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])</span>
  468 + <span class="comment">%set(fig_hdl, 'PaperSize', [fig_w fig_h]);</span>
  469 + set(fig_hdl, <span class="string">'PaperPositionMode'</span>, <span class="string">'auto'</span>);
  470 + print(fig_hdl,<span class="string">'-dpng'</span>,fname)
  471 + <span class="comment">% saveas(fig_hdl,fname,'png')</span>
  472 + <span class="comment">%crop([fname '.png'])</span>
  473 +<span class="keyword">end</span>
  474 +
  475 +
  476 +
  477 +<span class="keyword">for</span> i = 2:12
  478 + fig_hdl = figure(i);
  479 + <span class="keyword">if</span> save_flag
  480 + fname = [<span class="string">'Phases_'</span> site <span class="string">'_pair'</span> num2str(i-1) <span class="string">'-'</span> num2str(i)];
  481 +<span class="comment">% fig_w = 8;%cm</span>
  482 +<span class="comment">% fig_h = 5;%cm</span>
  483 +<span class="comment">% set(fig_hdl,'PaperUnits','centimeters')</span>
  484 +<span class="comment">% set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])</span>
  485 +<span class="comment">% set(fig_hdl, 'PaperSize', [fig_w fig_h]);</span>
  486 +<span class="comment">% set(fig_hdl, 'PaperPositionMode', 'auto');</span>
  487 +
  488 + print(fig_hdl,<span class="string">'-dpng'</span>,fname)
  489 + <span class="comment">% crop([fname '.png'])</span>
  490 + <span class="keyword">end</span>
  491 +<span class="keyword">end</span>
  492 +
  493 +
  494 +
  495 +fig_hdl = figure(13);
  496 +<span class="keyword">if</span> save_flag
  497 + fname = [<span class="string">'PhaseDifference_'</span> site <span class="string">'_calibrations'</span>];
  498 +<span class="comment">% fig_w = 8;%cm</span>
  499 +<span class="comment">% fig_h = 5;%cm</span>
  500 +<span class="comment">% set(fig_hdl,'PaperUnits','centimeters')</span>
  501 +<span class="comment">% set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])</span>
  502 +<span class="comment">% set(fig_hdl, 'PaperSize', [fig_w fig_h]);</span>
  503 + print(fig_hdl,<span class="string">'-dpng'</span>,fname)
  504 + <span class="comment">% crop([fname '.png'])</span>
  505 +<span class="keyword">end</span>
  506 +
  507 +fig_hdl = figure(14);
  508 +<span class="keyword">if</span> save_flag
  509 + fname = [<span class="string">'PhaseCalibration_'</span> site];
  510 +<span class="comment">% fig_w = 8;%cm</span>
  511 +<span class="comment">% fig_h = 5;%cm</span>
  512 +<span class="comment">% set(fig_hdl,'PaperUnits','centimeters')</span>
  513 +<span class="comment">% set(fig_hdl,'PaperPosition',[0 0 fig_w fig_h])</span>
  514 +<span class="comment">% set(fig_hdl, 'PaperSize', [fig_w fig_h]);</span>
  515 + print(fig_hdl,<span class="string">'-dpng'</span>,fname)
  516 + <span class="comment">% crop([fname '.png'])</span>
  517 +<span class="keyword">end</span>
  518 +</pre><pre class="codeoutput">
  519 +h6 =
  520 +
  521 + 3.3770e+03
  522 +
  523 +
  524 +h7 =
  525 +
  526 + 3.3780e+03
  527 +
  528 +</pre><p class="footer"><br>
  529 + Published with MATLAB&reg; 7.14<br></p></div><!--
  530 +##### SOURCE BEGIN #####
  531 +clear all
  532 +close all
  533 +
  534 +% calibrations made with the Helzel system
  535 +
  536 +
  537 +%cd('/home/bandma01/Documents/WERAs/')
  538 +cd('/home/bandma01/hfr-data-processing/calibrations/')
  539 +% datadir = '/run/media/bandma01/SeagateWireless/MEOPAR/DATA/HFRs/WERAs/';
  540 +% cd(datadir)
  541 +
  542 +% site = 'PAO';
  543 +site = 'PAB';
  544 +
  545 +cd(site)
  546 +
  547 +% calibrations: WALK or BOAT1 or BOAT2
  548 +for sys = 1:3
  549 +
  550 + if sys == 1
  551 + systm = 'WALK';
  552 + elseif sys == 2
  553 + systm = 'BOAT1';
  554 + elseif sys == 3
  555 + systm = 'BOAT2';
  556 + end
  557 +
  558 + clear J I E Phase calib_file calib_hr calib_jd calib_min calib_time calib_value
  559 + clear calib_year d data delta_phase_obs delta_phase_th gps_data index_timecalib
  560 + clear theta_pairs
  561 +
  562 + if strcmp(site,'PAO') == 1
  563 +
  564 + % load('HelzelPAO_Track_29-OCT-14 183138.mat')
  565 +
  566 + elseif strcmp(site,'PAB') == 1
  567 +
  568 + if strcmp(systm,'WALK') == 1
  569 + load('Track_05-JUL-16 150101.mat')
  570 + elseif strcmp(systm,'BOAT1') == 1
  571 + load('FirstCircle_BOAT_Track_05-JUL-16.mat')
  572 + elseif strcmp(systm,'BOAT2') == 1
  573 + load('SecondCircle_BOAT_Track_05-JUL-16.mat')
  574 + end
  575 +
  576 + end
  577 +
  578 +
  579 + % phases and time of calibration
  580 + if strcmp(site,'PAO') == 1
  581 + %cd('2014302')
  582 + elseif strcmp(site,'PAB') == 1
  583 + cd('2016187')
  584 + end
  585 +
  586 + % list files
  587 + clear d
  588 + d = dir('*.mat'); % returns a structure with name, date, bytes, isdir, datenum
  589 + Nfiles = length(d); % nombre de fichiers a analyser
  590 +
  591 +
  592 + % look at the names of files and extract information
  593 + for c = 1:Nfiles
  594 + calib_year(c) = str2double(d(c).name(1:4));
  595 + calib_jd(c) = str2double(d(c).name(5:7));
  596 + calib_hr(c) = str2double(d(c).name(8:9));
  597 + calib_min(c) = str2double(d(c).name(10:11));
  598 +
  599 + load(d(c).name,'phase')
  600 +
  601 + Phase(:,c) = phase;
  602 + end
  603 + clear phase
  604 +
  605 + % for each time stamp (column) compute phase difference between 2
  606 + % consecutive antennas (row)
  607 + delta_phase_obs = modsym(diff(Phase,1),180);
  608 +
  609 + To = datenum(2016,01,01,0,0,0);
  610 +
  611 +% % add 20 sec to center calibration time
  612 +% calib_time = To + calib_jd -1 + calib_hr./24 + calib_min/24/60 + 20/24/3600;
  613 +% % datevec(calib_time)
  614 +%
  615 +% clear calib_year calib_jd calib_hr calib_min
  616 +
  617 + calib_time = To + calib_jd -1 + calib_hr./24 + calib_min/24/60;
  618 +
  619 + if strcmp(site,'PAO') == 1
  620 + % calib time 2min34sec late wrt GPS
  621 + calib_time = calib_time + 2/24/60 + 34/24/3600;
  622 + end
  623 + disp(['the acquisition of calibration started at ' datestr(calib_time(1)) ' and ended up at ' datestr(calib_time(end))])
  624 +
  625 + %% separate between calibrations made walking and with the boat
  626 + % date of new calibration files
  627 + % walk 2016-07-05 13:53 and 14:58 UTC
  628 + % boat 1st circle 2016-07-05 18:20 and 19:18 UTC
  629 + % boat 2nd circle 2016-07-05 19:20 and 20:11 UTC
  630 +
  631 + if strcmp(systm,'WALK') == 1
  632 + T1 = datenum(2016,07,05,13,53,00);
  633 + T2 = datenum(2016,07,05,14,58,00);
  634 + elseif strcmp(systm,'BOAT1') == 1
  635 + T1 = datenum(2016,07,05,18,20,00);
  636 + T2 = datenum(2016,07,05,19,18,00);
  637 + elseif strcmp(systm,'BOAT2') == 1
  638 + T1 = datenum(2016,07,05,19,20,00);
  639 + T2 = datenum(2016,07,05,20,11,00);
  640 + end
  641 +
  642 + clear J K
  643 + J = find(calib_time < T1);
  644 + K = find(calib_time > T2);
  645 + J = [J K]; clear K
  646 +
  647 + % remove unwanted data to avoid confusion!
  648 +% calib_year(J)=[];
  649 +% calib_jd(J)=[];
  650 +% calib_hr(J)=[];
  651 +% calib_min(J)=[];
  652 + Phase(:,J)=[];
  653 + delta_phase_obs(:,J)=[];
  654 + calib_time(J)=[];
  655 + d(J)=[];
  656 + clear c Nfiles
  657 +
  658 + %+++++++++++++++++++++++++++++++++++++++++
  659 + % for each calibration file (series# 1)
  660 + %+++++++++++++++++++++++++++++++++++++++++
  661 + for calib_file = 1:length(calib_time)
  662 +
  663 + time_calib = calib_time(calib_file);
  664 + % datevec(time_calib)
  665 +
  666 + % find corresponding file in gps data
  667 + clear I
  668 + I = find(time_calib <= data(:,1) );
  669 + % take closest point for now
  670 + I = I(1);
  671 +
  672 + % datevec(data(I,1))
  673 + % datevec(time_calib)
  674 +
  675 + index_timecalib(calib_file) = I;
  676 + end
  677 +
  678 + % check
  679 + %datevec(data(index_timecalib,1))
  680 +
  681 + % isolate gps calibration data
  682 + gps_data = data(index_timecalib,:); % time lat lon
  683 +
  684 + figure(1)
  685 + hold on
  686 + if sys == 1
  687 + h1(sys) = plot(gps_data(:,2),gps_data(:,3),'.-','displayName','gps WALK');
  688 + elseif sys == 2
  689 + h1(sys) = plot(gps_data(:,2),gps_data(:,3),'.-','displayName','gps BOAT circle 1');
  690 + elseif sys == 3
  691 + h1(sys) = plot(gps_data(:,2),gps_data(:,3),'.-','displayName','gps BOAT circle 2');
  692 + end
  693 +
  694 + grid on
  695 +
  696 + if strcmp(systm,'WALK') == 1
  697 + set(h1(1),'marker','x','color','b')
  698 + elseif strcmp(systm,'BOAT1') == 1
  699 + set(h1(2),'marker','.','color',[0 0.5 0])
  700 + elseif strcmp(systm,'BOAT2') == 1
  701 + set(h1(3),'marker','*','color',[0.6 0.2 0])
  702 + end
  703 +
  704 + % load antenna coordinates
  705 + %
  706 + %load([site '_Rx_coords.mat'])
  707 +
  708 +
  709 + if strcmp(site,'PAO') == 1
  710 + cd('/home/bandma01/hfr-data-processing/calibrations/')
  711 + load(['Survey_2014_' site])
  712 + lon = PAO.lon;
  713 + lat = PAO.lat;
  714 + elseif strcmp(site,'PAB') == 1
  715 + cd ..
  716 + load([site '_Rx_coords.mat'])
  717 + end
  718 +
  719 + hold on
  720 + h2b(sys) = plot(lon,lat,'dk','displayName','antennas positions');
  721 + grid on
  722 + if strcmp(systm,'WALK') == 1
  723 + %legend('gps Walk','antennas positions','location','best')
  724 + legend([h2b h1],'location','best')
  725 + elseif strcmp(systm,'BOAT1') == 1
  726 + legend([h2b(1) h1],'location','best')
  727 + % legend('gps Walk','antennas positions','gps Boat circle 1','antennas positions','location','best')
  728 + elseif strcmp(systm,'BOAT2') == 1
  729 + legend([h2b(1) h1],'location','best')
  730 + % legend('gps Walk','antennas positions','gps Boat circle 1','antennas positions','gps boat circle 2','antennas positions','location','best')
  731 + end
  732 + xlabel('longitude')
  733 + ylabel('latitude')
  734 + title(['site = ' site ])
  735 +
  736 +% continue
  737 +
  738 + %% compute theoretical phase values
  739 + % wavelength
  740 + lambda = 3E8/16.15E6;
  741 +
  742 +
  743 + % initialize variables
  744 + theta_pairs = nan(11,length(gps_data));
  745 +
  746 + % for each gps point
  747 + for pt = 1:length(gps_data)
  748 +
  749 + % compute distance between GPS and antennas
  750 + for ant = 1:12
  751 + clear E W
  752 + [E,W] = lonlat2km(gps_data(pt,2),gps_data(pt,3),lon(ant),lat(ant));
  753 + ddist = sqrt(E.^2+W.^2);
  754 + r(ant) = ddist*1000; %meters
  755 + end
  756 +
  757 +
  758 + % for each pair of antennas
  759 + for p = 1:11
  760 +
  761 + %delta_r
  762 + delta_r(p) = r(p+1) - r(p);
  763 +
  764 + % theoretical phase
  765 + delta_phase_th(p,pt) = modsym((delta_r(p)*360./lambda),180);
  766 +
  767 +
  768 + % position of middle of pairs
  769 + clear lon_mid lat_mid
  770 + lon_mid = (lon(p+1)+lon(p))/2;
  771 + lat_mid = (lat(p+1)+lat(p))/2;
  772 +
  773 + % distance between middle of pair and one of the antenna in the pair
  774 + % clear E W
  775 + % [E,W] = lonlat2km(lon_mid,lat_mid,lon(p),lat(p));
  776 + % ddist = sqrt(E.^2+W.^2);
  777 + % a(p) = ddist*1000; %meters
  778 +
  779 + % distance between middle of pair and GPS pt
  780 + clear ddsit phasAng
  781 + [ddist,phasAng] = sw_dist([gps_data(pt,3) lat_mid],[gps_data(pt,2) lon_mid],'km');
  782 + r_mid(p) = ddist*1000; %meters
  783 + zeta(p) = phasAng;
  784 +
  785 +
  786 + % distance between pair of antennas
  787 + clear ddsit phasAng
  788 + [ddist,phasAng] = sw_dist([lat(p) lat(p+1)],[lon(p) lon(p+1)],'km');
  789 + phi(p) = phasAng;
  790 +
  791 +
  792 + % azimuth
  793 + theta_pairs(p,pt) = modsym(zeta(p)-phi(p)+90,180);
  794 + %rad2deg(theta_pairs)
  795 + end
  796 + end
  797 +
  798 +
  799 +
  800 + % for each pair of antenna
  801 + for p = 1:11
  802 +
  803 + figure(p+1)
  804 +
  805 + h2(sys,p) = plot(theta_pairs(p,:),delta_phase_th(p,:),'-','displayName','theoretical');
  806 + hold on
  807 + if sys == 1
  808 + h3(sys,p) = plot(theta_pairs(p,:),-delta_phase_obs(p,:),'o','displayName','WALK calib');
  809 + elseif sys == 2
  810 + h3(sys,p) = plot(theta_pairs(p,:),-delta_phase_obs(p,:),'o','displayName','BOAT circle 1 calib');
  811 + elseif sys == 3
  812 + h3(sys,p) = plot(theta_pairs(p,:),-delta_phase_obs(p,:),'o','displayName','BOAT circle 2 calib');
  813 + end
  814 + grid on
  815 + xlabel('incidence (degrees)')
  816 + ylabel('phase difference (degrees)')
  817 +
  818 + if strcmp(systm,'WALK') == 1
  819 + set(h2(sys,p),'color','b')
  820 + set(h3(sys,p),'color','r','marker','o')
  821 + legend([h2(1,p) h3(1,p)],'location','best')
  822 + %legend('theoretical','Walk','location','best')
  823 + elseif strcmp(systm,'BOAT1') == 1
  824 + set(h2(sys,p),'color','b')
  825 + set(h3(sys,p),'color','r','marker','sq','markerfacecolor','r')
  826 + %legend('theoretical','Boat circle 1','location','best')
  827 + legend([h2(1,p) h3(1,p) h3(2,p)],'location','best')
  828 + elseif strcmp(systm,'BOAT2') == 1
  829 + set(h2(sys,p),'color','b')
  830 + set(h3(sys,p),'color',[0 0.5 0],'marker','sq','markerfacecolor',[0 0.5 0])
  831 + %legend('theoretical','Boat circle 2','location','best')
  832 + legend([h2(1,p) h3(1,p) h3(2,p) h3(3,p)],'location','best')
  833 + end
  834 + title([site ' - antennas ' num2str(p) '-' num2str(p+1)])
  835 +
  836 +
  837 + % compute calibration value for antennas pairs
  838 + % calib_value(p) = nanmean(modsym(-delta_phase_obs(p,:)-delta_phase_th(p,:),180),2);
  839 + calib_value(p) = nanmean(modsym(delta_phase_th(p,:)+delta_phase_obs(p,:),180),2);
  840 +
  841 +
  842 + % compute calibration value for each antenna
  843 + phi_ant(p+1,sys) = sum(calib_value);
  844 +
  845 + end
  846 +
  847 + %continue
  848 +
  849 + P = ones(16,3);
  850 + P(13:16,:) = 0;
  851 + P(1:12,3) = phi_ant(:,1);
  852 +
  853 +
  854 + % save file (we only keep results from Helzel
  855 +% datasav = [datadir site];
  856 +% cd(datasav)
  857 +% if strcmp(site,'PAO') == 1
  858 +% dlmwrite('calibration_pao_29OCT2014.wera',P,'delimiter','\t','precision','%.7f')
  859 +% elseif strcmp(site,'PAB') == 1
  860 +% dlmwrite('calibration_pab_30OCT2014.wera',P,'delimiter','\t','precision','%.7f')
  861 +% end
  862 +
  863 +
  864 + figure(13)
  865 + hold on
  866 + if sys == 1
  867 + h4(sys) = plot(calib_value,'.-','displayname','WALK','color','r','marker','o');
  868 + elseif sys == 2
  869 + h4(sys) = plot(calib_value,'.-','displayname','BOAT circle 1','color','k','marker','sq');
  870 + elseif sys == 3
  871 + h4(sys) = plot(calib_value,'.-','displayname','BOAT circle 2','color',[0 0.5 0],'marker','d');
  872 + end
  873 + legend(h4)
  874 + grid on
  875 + xlabel('pairs of antennas')
  876 + ylabel(['averaged phase differences between theoretical',10 ,'and observed values (degrees)'])
  877 + title(site)
  878 +
  879 +
  880 +
  881 + figure(14)
  882 + hold on
  883 + if sys == 1
  884 + h5(sys) = plot([1:12],phi_ant(:,sys),'displayname','WALK','color','r','marker','o');
  885 + elseif sys == 2
  886 + h5(sys) = plot([1:12],phi_ant(:,sys),'displayname','BOAT circle 1','color','k','marker','sq');
  887 + elseif sys == 3
  888 + h5(sys) = plot([1:12],phi_ant(:,sys),'displayname','BOAT circle 2','color',[0 0.5 0],'marker','d');
  889 + end
  890 + grid on
  891 + legend(h5)
  892 + xlabel('antenna number')
  893 + ylabel('calibration phase to be applied (degrees)')
  894 + title(site)
  895 +
  896 + if strcmp(site,'PAO') == 1
  897 + break
  898 + end
  899 +
  900 +end
  901 +
  902 +
  903 +% add current calibration phases
  904 +if strcmp(site,'PAO') == 1
  905 + load calibration_pao.wera
  906 + figure(14)
  907 + hold on
  908 + plot([1:12],calibration_pao(1:12,3),'b.-')
  909 + legend('Helzel','original calibration','location','northwest')
  910 +elseif strcmp(site,'PAB') == 1
  911 + load calibration_pab.wera
  912 + figure(14)
  913 + hold on
  914 + h6 = plot([1:12],calibration_pab(1:12,3),'b.-','displayname','primary calibrations')
  915 +
  916 + load calibration_pab_30OCT2014.wera
  917 + figure(14)
  918 + hold on
  919 + h7 = plot([1:12],calibration_pab_30OCT2014(1:12,3),'m.-','displayname','calibrations 30OCT2014')
  920 + legend([h5 h6 h7],'location','best')
  921 +
  922 +end
  923 +
  924 +
  925 +
  926 +
  927 +
  928 +
  929 +% save figures
  930 +save_flag = 0;