%%%  plan_deploy13.m  %%%%%
%
% Used to generate Buoy Tasking Orders   BTO
%
%  Rick Lumpkin & Eric D'Asaro, August 2008

clear all;

[filename, pathname] = uigetfile('ARRAY*.m', 'Select ARRAY file'); %,'multiselect','off');
arrayfile=[pathname filename];
run(arrayfile)


% figure out number of lines in drifter array
Nlines=1;
for i=2:50;
    if exist(['line',num2str(i)]) Nlines=Nlines+1; end;
end;

t_dep=juliand(dep);


% load hurricane locations and wind fields
% (downloaded from NHC web site)
loadhurrwinds
clf;

% load image of absolute sea level from Navy NLOM SSH nowcast
% for Inter-American Seas 
%   http://www7320.nrlssc.navy.mil/global_nlom32/navo/IAS.gif
if exist('SSHplot')
    if SSHplot==1
        %eval(['!wget -N -r -nd --retr-symlinks ' ...
        %    '"http://www7320.nrlssc.navy.mil/global_nlom32/navo/IAS.gif"']); 
        [ssh,colmap]=imread('IAS.gif','gif');
        ssh=flipud(ssh(55:530,31:730));
        lonssh=-100:(-55+100)/699:-55;
        latssh=5:(31-5)/475:31;
        imagesc(lonssh,latssh,ssh);
        colormap(colmap);
        set(gca,'ydir','normal');
        hold on;
    end;
end;

% Plot coastline and bathymetry
%coasts(3,1,[],'k');
hold on;
load bathycarib;
if SSHplot~=1  % No SSH
    [c,h]=contour(lon2,lat2,z2,[200 200],'k');
    clabel(c,h,'FontSize',8);
    [c,h]=contour(lon2,lat2,z2,[1500 1500],':k');
    clabel(c,h,'FontSize',8);
else   % if SSH contour another way
    c=contourc(lon2,lat2,z2,[200 1500]);
    ct=c;
    while size(ct,2)>0
        ns=ct(2,1);
        level=ct(1,1);
        h=plot(ct(1,2:ns+1),ct(2,2:ns+1),'k-','Color',[1 1 1]*0.8);hold on;
        ct=ct(:,ns+2:end);
    end
end

%EEZplot;

set(gca,'ydir','normal')
hold on;

% interpolate center location and wind radii for deployment time
lon_hurr_dep=interp1(t,lon,t_dep,'linear');
lat_hurr_dep=interp1(t,lat,t_dep,'linear');
for j=1:4;
    r64_dep(1,j)=interp1(t,r64(:,j),t_dep,'linear');
    r50_dep(1,j)=interp1(t,r50(:,j),t_dep,'linear');
    r34_dep(1,j)=interp1(t,r34(:,j),t_dep,'linear');
end;
plot(lon,lat,'b','linewidth',2);
hold on;
windpatch(t_dep,lon_hurr_dep,lat_hurr_dep,r64_dep,r50_dep,r34_dep);
axis([lon_hurr_dep-6*max(max(r34))/60 lon_hurr_dep+2*max(max(r34)/60) ...
       lat_hurr_dep-4*max(max(r34))/60 lat_hurr_dep+4*max(max(r34))/60]);
hold on;
plot(lon_hurr_dep,lat_hurr_dep,'b+','MarkerSize',12,'LineWidth',2);

% plot center and gale-force (34kt) wind radius for deployment time + 24h
lond2=interp1(t,lon,t_dep+1,'linear');
latd2=interp1(t,lat,t_dep+1,'linear');
plot(lond2,latd2,'xk')
if (length(find(isfinite(lond2))))
    for j=1:4;
        r34d2(1,j)=interp1(t,r34(:,j),t_dep+.1,'linear');
    end;
    theta=0:pi/16:2*pi;
    R=NaN*ones(size(theta));
    del=abs(theta-pi/4); % NE
    q=find(del==min(del));    R(q)=r34d2(1);
    del=abs(theta-7*pi/4);   % SE
    q=find(del==min(del));   R(q)=r34d2(2);    
    del=abs(theta-5*pi/4);   % SW
    q=find(del==min(del));   R(q)=r34d2(3);
    del=abs(theta-3*pi/4);  % NW
    q=find(del==min(del));  R(q)=r34d2(4);
    q=find(theta==0 | theta==max(theta));
    R(q)=( r34d2(1)+r34d2(2) )/2;
    q=find(isnan(R));
    q2=find(isfinite(R));  
    R(q)=interp1(theta(q2),R(q2),theta(q),'linear');
    x=R.*cos(theta); x=[x x(1)];
    y=R.*sin(theta); y=[y y(1)];
    plot(lond2+x/60,latd2+y/60,'--k')
end
hold off;

% plot all storm predictions
if Onestorm~=1
    hold on
    windpatch(t,lon,lat,r64,r50,r34);
    hold on
end

if (~minor_adjust)
    Nlines=input('Number of lines in array? ');
end;

% store the Lon, Lat, Line Number, Course, Sep. Distance,
% and buoy type for each drop point.  "Type" is 10 values
% across, allowing multiple types of instruments at each
% drop point (e.g., a drifter and two floats), with NaN
% as the placeholder.
Lon=[]; Lat=[]; Line=[]; Course=[]; Dist=[]; Type=[];
for i=1:Nlines


    eval(['linefile=line',num2str(i),';']); 

    % place line of drifters closest to storm
    if (i==1)
        if (~minor_adjust)
            disp('click on center of the line closest to storm,')
            disp('  and again to define direction of the line (and of the flight)')
            [lon_dep,lat_dep]=ginput(2)
        else
            lon_dep=[linefile(1).lon(1) linefile(1).lon(2)];
            lat_dep=[linefile(1).lat(1) linefile(1).lat(2)]; 
        end;
    else
        %recip=input('Is the next line on the reciprocal course (y/n)','s');
        %if (strcmp(recip,'y'))

        if (~minor_adjust)
            fprintf('click on center of line %d,\n',i)
            disp('  and again to define direction of line/flight')
            [lon_dep,lat_dep]=ginput(2)
        else 
            lon_dep=[linefile(1).lon(1) linefile(1).lon(2)];
            lat_dep=[linefile(1).lat(1) linefile(1).lat(2)]; 
        end;
    end;

    [dist,ang]=sw_dist(lat_dep,lon_dep,'nm');
    ang=round(ang);
    % convert to degrees true, compass convention.  sw_dist angle is
    % positive counter-clockwise, with 0 being due east.  compass
    % convention is positive clockwise, with 0 being due north
    theta=-(ang-90);  if (theta<0) theta=theta+360; end; 

    if (~minor_adjust)
        Ndrift=input('How many drifters in this line? ');
        d=input('Separation distance (nm) (default=30) ');
        if (~length(d)) d=30; end;
    else
        Ndrift=length(linefile);
        d=linefile(1).spacing;
    end;


    % the line of drifters will be (Ndrift-1)*d nautical miles long.
    % starting point is at recripocal ang+180

    Lon=[Lon; lon_dep(1)- (Ndrift-1)/2*d*(diff(lon_dep)/dist)];
    Lat=[Lat; lat_dep(1)- (Ndrift-1)/2*d*(diff(lat_dep)/dist)];
    Line=[Line; i];
    Course=[Course; theta];
    Dist=[Dist;d];
    type=NaN*ones(1,10);
    for n=1:length(linefile(1).buoytype)
        type(n)=linefile(1).buoytype(n);
    end;
    Type=[Type;type];
    for j=2:Ndrift;
        Lon=[Lon; Lon(length(Lon))+d*(diff(lon_dep)/dist)];
        Lat=[Lat; Lat(length(Lat))+d*(diff(lat_dep)/dist)];
        Line=[Line; i];
        Course=[Course; theta];
        Dist=[Dist;d];
        type=NaN*ones(1,10);
        for n=1:length(linefile(j).buoytype)
            type(n)=linefile(j).buoytype(n);
        end;
        Type=[Type;type];
    end;
    hold on;
    plot(Lon,Lat,'ok');
    hold off;

end;    % looping through drop points, storing values


hold on;
% plot estimated storm location at 12h increments before and after t_dep
tt=t_dep-0.5:0.5:t_dep+2;
q=find(tt<min(t) | tt>max(t) | tt==t_dep); tt(q)=[];
llat=interp1(t,lat,tt,'linear');
llon=interp1(t,lon,tt,'linear');
plot(llon,llat,'.k','markersize',20);
q=find(tt>t_dep);
for i=1:length(q);
    txt=['+',num2str(12*i),'h'];
    text(llon(q(i)),llat(q(i)),txt,'VerticalAlignment','top','FontSize',9,'Color','w',...
        'FontWeight','Bold');
    text(llon(q(i)),llat(q(i)),txt,'VerticalAlignment','top','FontSize',9,'Color','b');

end;

% plot planned array
for i=1:Nlines;
    q=find(Line==i);
    plot(Lon(q),Lat(q),'k');
end;
for i=1:length(Lon);
    if (length(find(Type(i,:)==0)))  % MINIMET
        plot(Lon(i),Lat(i),'ok','markersize',12,'markerfacecolor',[1 0 1]);
    end;
    if (length(find(Type(i,:)==1)))    % ADOS
        plot(Lon(i),Lat(i),'ok','markersize',12,'markerfacecolor','w');
    end;
    if (length(find(Type(i,:)==2)))    % gas float
        offset=(Course(i)-90)*2*pi/360;
        text(Lon(i)+Dist(i)/2/60*sin(offset), ...
            Lat(i)+Dist(i)/2/60*cos(offset),'g','fontsize',12);
    end;
    if (length(find(Type(i,:)==3)))    % EM float
        offset=(Course(i)+90)*2*pi/360;
        text(Lon(i)+Dist(i)/2/60*sin(offset), ...
            Lat(i)+Dist(i)/2/60*cos(offset),'e','fontsize',12);
    end;    
    text(Lon(i),Lat(i),num2str(i),'horizontalalignment','center');
end;
hold off;


flight_dist=sw_dist([takeoff(1) Lat(1)], ...
    [takeoff(2) Lon(1)],'nm');
% hours between plane takeoff and first deployment 
flight_hours=.25*round(4*flight_dist/AirSpeed);   


% create tasking order, assuming UTC is current time plus 4 hours
fprintf(1,'Assuming UTC is current time plus %4.1f hours\n',LOCAL2GMT)
time_now=datestr(now+LOCAL2GMT/24);
gd_dep=datestr(t_dep + 722815);
gd_depP=datestr(t_dep+4/24 + 722815);
gd_depTO=datestr(t_dep-flight_hours/24 + 722815);

Task={'BUOY TASKING ORDER'};
Task(2)={['DTG: ',time_now(1:2),'/',time_now(13:14), ...
        time_now(16:17),' ',time_now(4:6),' ',time_now(10:11)]};
Task(3)={['VALID FOR DEPLOYMENTS BETWEEN ',gd_dep(1:2),'/',gd_dep(13:14), ...
        gd_dep(16:17),'Z ',gd_dep(4:6),' ',gd_dep(10:11),' AND ', ...
        gd_depP(1:2),'/',gd_depP(13:14), ...
        gd_depP(16:17),'Z ',gd_depP(4:6),' ',gd_depP(10:11)]};
Task(4)={['TAKEOFF TIME ',gd_depTO(1:2),'/',gd_depTO(13:14), ...
        gd_depTO(16:17),'Z ',gd_depTO(4:6),' ',gd_depTO(10:11)]};
Task(5)={' '};
counter=1;
% loop through lines (index i)
for i=1:Nlines

    % construct header line.  Example:
    % LINE 1: SP 26 45N 92 44W  EP 23 15N 92 44W HEADING 180T SPACING/RADIUS 30/5nm
    q=find(Line==i);
    deg=floor(Lat(q(1))); 
    mn=round(60*(Lat(q(1))-floor(Lat(q(1)))));
    pos_txt=['SP ',num2str(deg),' '];
    if (mn<10) pos_txt=[pos_txt,'0']; end;
    pos_txt=[pos_txt,num2str(mn),'N'];
    deg=floor(abs(Lon(q(1))));
    mn=round(60*(abs(Lon(q(1)))-floor(abs(Lon(q(1))))));
    pos_txt=[pos_txt,' ',num2str(deg),' '];
    if (mn<10) pos_txt=[pos_txt,'0']; end;
    pos_txt=[pos_txt,num2str(mn),'W  '];

    deg=floor(Lat(q(length(q)))); 
    mn=round(60*(Lat(q(length(q)))-floor(Lat(q(length(q))))));
    pos_txt=[pos_txt,'EP ',num2str(deg),' '];
    if (mn<10) pos_txt=[pos_txt,'0']; end;
    pos_txt=[pos_txt,num2str(mn),'N'];
    deg=floor(abs(Lon(q(length(q)))));
    mn=round(60*(abs(Lon(q(length(q))))-floor(abs(Lon(q(length(q)))))));
    pos_txt=[pos_txt,' ',num2str(deg),' '];
    if (mn<10) pos_txt=[pos_txt,'0']; end;
    pos_txt=[pos_txt,num2str(mn),'W'];

    Task(length(Task)+1)={['LINE ',num2str(i),': ',pos_txt, ...
                ' HEADING ',num2str(Course(q(1))),'T  SPACING/RADIUS ', ...
                num2str(Dist(q(1))),'/5nm']};
    %  end, constructing header line

    %  individual lines of BTO for the instruments.  Example:
    %  "1:  MINIMET"  or "3:  ADOS, Lagrangian 66" 
    for j=1:length(q);

        buoyline=[num2str(counter),': '];
        % loop through all 10 "Type" entries, filling in text
        %  for non-NaN values.
        for k=1:10;

            if (k>1 & isfinite(Type(counter,k)))
                buoyline=[buoyline,', '];
            end;

            if (Type(counter,k)==0)
                buoyline=[buoyline,['MINIMET-',num2str(BoxNumM(1))]];
                BoxNumM(1)=[];
            elseif (Type(counter,k)==1)
                buoyline=[buoyline,['ADOS-',num2str(BoxNumA(1))]];
                BoxNumA(1)=[];
            elseif (Type(counter,k)==2) % Gas float
                buoyline=[buoyline,['L-',num2str(BoxNumGas(1))]];
                BoxNumGas(1)=[];
            elseif (Type(counter,k)==3) % EM float
                buoyline=[buoyline,['L-',num2str(BoxNumEM(1))]];
                BoxNumEM(1)=[];
            end;

        end;
        Task(length(Task)+1)={buoyline};
        counter=counter+1;
    end;    % loop over instrument(s) line "1: MINIMET", "2: ADOS, Lagrangian 66"
    Task(length(Task)+1)={' '};
end;
Task(length(Task)+1)={['END BUOY TASKING ORDER DTG: ', ...
            time_now(1:2),'/',time_now(13:14), ...
            time_now(16:17),' ',time_now(4:6),' ',time_now(10:11)]};


fprintf('\n\n\n');
for i=1:length(Task);
    disp(char(Task(i)))
end;

% FIX UP PLOT
minlon=min(LonP); maxlon=max(LonP);
minlat=min(LatP); maxlat=max(LatP);
axis([minlon maxlon minlat maxlat]);
latdist=1;
londist=cos((maxlat+minlat)/2*2*pi/360);
set(gca,'dataaspectratio', [latdist londist 1]);


% mapaxis seems screwed up
% mapaxis([min(LonP) max(LonP) min(LatP) max(LatP)])

text(minlon+.1,minlat+.3,'White: ADOS (t-chain).  Magenta: MIMIMET');
gd_pred=datestr(min(t)+ 722815);
text(max(LonP)-.1,max(LatP)-.3,[gd_pred(1:6),' ', ...
         gd_pred(13:14),gd_pred(16:17),'Z forecast'], ...
     'horizontalalignment','right');
hold on;
xlabel('Degrees East Longitude');
ylabel('Degrees North Latitude');

if exist('gus_lon');
    hold on;
    ID=unique(gus_id);
    for i=1:length(ID);
        q=find(gus_id==ID(i));
        plot(gus_lon(q),gus_lat(q),'k');
        % plot most recent fix, if less than 1 day old
        qq=find(gus_t(q)==max(gus_t(q)));
        if (gus_t(q(qq))>max(gus_t)-1) 
            plot(gus_lon(q(qq)),gus_lat(q(qq)), ...
                'ow','markersize',6,'markerfacecolor','b')
        end;
    end;
end;
hold off;

%%%%%%%%%%%%% END plan_deply  %%%%%%%%%%%%
