% routine to post coherence, phase, & gain on a map bewteen sites
% for various frequency levels
% (Note: Input is obtained by running wind/mfiles/cspec.m.)

% loop for nf frequencies ranges
freq=[0 2;2 10];
[nf,nc]=size(freq);
for j=1:nf;
  
  %subplot by frequency
  subplot(nf,1,j);

  % make basemap
  georbat(1,0,0,2,1);
  axis([-70.5 -66 40 42]);
  lllabel(60,60,40,42,-70,-66);
  
  % load lat & lon
  load /home3/ocn/jmanning/wind/stations.dat
  s=stations;
  format bank

  % how many pairs to plot
  hmp=3;

  % loop through pairs
  for k=1:hmp;

    case=input('Enter case: ','s');

    % load coh&ph info
    eval(['load /home3/ocn/jmanning/wind/mfiles/cspecdat/cspec',case,'.dat']);
    eval(['c=cspec',case,';']);


    % what pair is this?
    s
    wp=input('What pair is this?');

    % plot the two points
    plot([s(wp(1),2) s(wp(2),2)],[s(wp(1),3) s(wp(2),3)],...
      'w.','markersize',60);
    for i=1:2;
      text(s(wp(i),2),s(wp(i),3),...
        num2str(s(wp(i),1)),...
        'horizontalalignment','center','color','k',...
        'fontsize',12);
    end;

    % plot line between points
    %plot([s(wp(1),2) s(wp(2),2)],[s(wp(1),3) s(wp(2),3)],...
    %  'color',['c',num2str(k)]);
    

    % calculate angle from pt one to two
    [distkm,bear]=dist(s(wp(1),3),s(wp(1),2),s(wp(2),3),s(wp(2),2));
     bear0=bear;
    if bear>180;
     [distkm,bear]=dist(s(wp(2),3),s(wp(2),2),s(wp(1),3),s(wp(1),2));
    end;

    % post relavent info after finding results in this freq range
    id=find(c(:,3)>=freq(j,1)&c(:,3)<=freq(j,2));
    text((s(wp(1),2)+s(wp(2),2))/2,(s(wp(1),3)+s(wp(2),3))/2,...
      [num2str(meannnan(c(id,7)),2),', ',...
       num2str(meannnan(c(id,5)),2),', ',...
       num2str(meannnan(c(id,4)),2)],...
      'horizontalalignment','center',...
      'verticalalignment','bottom',...
      'color',['c',num2str(k)],...
      'rotation',-1*(bear-90));

    % add an arrow for original direction
    [u,v]=sd2uv(5,bear0);
    arrowold((s(wp(1),2)+s(wp(2),2))/2,(s(wp(1),3)+s(wp(2),3))/2,...
       u,v,1,1,1,[-70.5 -66 40 42],['c',num2str(k)]);

    % title
    title(['Coherence, Phase, and Gain for ',num2str(freq(j,1)),'-',num2str(freq(j,2)),' day cycles'],...
        'fontsize',14);

    %fix axis
    if j==1;set(gca,'xticklabels',' ');xlabel(' ');
       set(gca,'position',[.1 .55 .8 .4]);
    else;
       set(gca,'position',[.1 .05 .8 .4]);
    end;


  end; % go get another freqeuncy; 
end;