function recip(label, froot, tr)
% function recip(label, froot, tr)
%
% Calculate calibration phase and amplitude
% from reciprocal ship tracks; we treat these as a closed
% loop and integrate around the circuit.
%
% estimates "goodness" of closed loop assumption by
%  1) computing sum of gap lengths over total length of circuit ( as %)
%  2) area enclosed by circuit over total length *** NOT IMPLEMENTED ***
%
% input needed:
%   *sm file from nav/ (cp or create symbolic link to this directory)
%
%   command line arguments:
%
%   label: a string containing title (usually cruise ID)
%
%   froot: a string with filename root of *.sm file (smoothr output)
%
%   tr is an array of size n X 2 containing the endpoints
%     of the timeranges which define the segments of the
%     loop, in decimal days; input form [t0 t1; t2 t3; ...]
%       FOR EXAMPLE, tr=[296.7742 297.4743; 322.4506 323.2527] for
%       two tracks which are "reciprocal" but are not continuous
%       in time; tr=[298.3648 299.2122] for a continous reciprocal
%
% output:
%   recip.out: text file containing  amplitude and phase values
%               and goodness of assumption estimation
% 96/03/08 JF
% 01/12/12 JF: changed fprintf statements to use FID
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

area_flag = 0;          % disable area calculation

sm = read_asc([froot,'.sm']);            % load *sm file

te = sm(:,1);                            % ensemble times
[nsm,nc]=size(te);

[ntr,nn]=size(tr);                      % ntr = # segments

% create a logical mask, the size of t, where value is set
%  to 1 for those times which fall within the timeranges
%  and have ship and doppler velocity values

seg_mask=zeros(nsm,1);                  % initialize use mask

for ii=1:ntr
   seg_mask=seg_mask | ((te>= tr(ii,1) & te<= tr(ii,2)) & (sm(:,2)<1000));
end
i_use=find(seg_mask);                   % indices to use
[ne,ns]=size(i_use);

% get complex U_ship and U_doppler for desired times
%  doppler (ud) is water vel relative to ship, so subtract us from
%  abs ref layer vel to get doppler

us=sm(i_use,4) + i*sm(i_use,5);
ud=sm(i_use,2) + i*sm(i_use,3) - us;

% get time differences; check for delta t > 0.0040 and
%  replace with median(t) for amp/phs calc

t=te(i_use);          % take only needed times;

dt=diff(t);
it=find(dt > 0.0040);
dt(it)=median(dt)*ones(size(it));

tq= [t; t(1)];    % need last-first in quality calc
dtq = diff(tq);
itq = find(dtq > 0.0040 | dtq < 0.0);  % first-last < 0!

% calculate phase
% first find integral over time of ud*conj(us)

iud_us=sum(ud(2:ne).*conj(us(2:ne)).*dt);

% use angle function to solve for phi, convert to degrees
%       atan restricts output to +/- pi/2

phi_rad = -atan(imag(iud_us)/real(iud_us));
phi = phi_rad*180/pi

% now calculate amp

btm=real(exp(i*phi_rad).*ud(2:ne).*conj(us(2:ne)).*dt);
top=real(-us(2:ne).*conj(us(2:ne)).*dt);
amp=sum(top)/sum(btm)

% get measure of quality of approximation:
%   percentage of loop where there are gaps
%
%
% add first position to end of array to "close" the loop
%  get length in meters between positons for total length

lat = (sm(i_use,7));
lat_loop = [lat; lat(1)];
[ni,nl] = size(lat_loop);
dlat = diff(lat_loop);
alat = lat_loop(2:ni) - dlat/2;
dy = lat_to_m(dlat, alat);
lon = (sm(i_use,6));
lon_loop = [lon; lon(1)];
dlon = diff(lon_loop);
alon = lon_loop(2:ni) - dlon/2;
dx = lon_to_m(dlon, alon);
ds = sqrt(dy.^2 + dx.^2);
L = sum(ds);                     % total length of loop

% now get sum of gap lengths
%  gaps identified in index array itq (includes last to
%   first position gap)

ds_gap = ds(itq);
G = sum(ds_gap);

% percent of loop that gaps occupy

pc_gap = (G/L) * 100

% print out amp and phase values and quality measures

Ts=tr(1,1);
Te=tr(ntr,2);

FID=fopen('recip.out','a');
fprintf(FID, [label '\n'] );
fprintf(FID, [' Time range %9.4f to %9.4f \n'], tr.');
fprintf(FID,  ['\n   Calculation done at ' pctime ]);
fprintf(FID, '\n   amp   = %7.4f ', amp);
fprintf(FID, '\n   phase = %7.4f', phi);
fprintf(FID, '\n   pcgap = %8.4f', pc_gap);
fprintf(FID,'\n');
fclose(FID);

FID=fopen('amp_ph.out','a');
fprintf(FID,'\n %7.4f   %7.4f', amp, phi);
fclose(FID);

