%Generate ampl. and phases from FES2012
%Inputs: (1) open.ll (generated by gen_fg.f90 with 1st 2 lines removed, 1 seg at a time);  
%        (2) *_FES2012_SLEV_nc3.nc (in this dir)
%Output: ap.out (for bctides.in; see below for order)
clear all; close all;
const={'S2','M2','N2','K2','K1','P1','O1','Q1'};
%const={'O1','K1','Q1','M2','S2','K2','N2'};
open=load('open.ll'); %ID,lon,lat of open bnd nodes
npt=size(open,1);

fid=fopen('ap.out','w');

miss_value=1.e10; %junk value
for i=1:length(const)
  ncid = netcdf.open([const{i} '_FES2012_SLEV_nc3.nc'],'NC_NOWRITE');
  vid=netcdf.inqVarID(ncid,'lat'); %input var./array name
  lat = netcdf.getVar(ncid, vid); %ascending order (-90,90)
  vid=netcdf.inqVarID(ncid,'lon');
  lon = netcdf.getVar(ncid, vid); %ascending order [0,360)
  nx=length(lon); ny=length(lat);

  vid=netcdf.inqVarID(ncid,'Ha'); %(nx,ny)
  amp0= netcdf.getVar(ncid, vid);
  amp0=amp0/100; %to meters
  vid=netcdf.inqVarID(ncid,'Hg');
  pha0= netcdf.getVar(ncid, vid); %degr
  %Put phases into [0,360)
  indx=find(pha0<0);
  pha0(indx)=pha0(indx)+360;
  indx=find((pha0<0 | pha0>360) & pha0<1.4);
  if(sum(indx) ~=0) 
    disp('Phase out of bound');
    pha0(indx)
    error('');
  end;

  ap=zeros(npt,2);
  fprintf(fid,'%s\n',const{i});
  for j=1:npt
    lon2=open(j,2); if(lon2<0); lon2=lon2+360; end;
    lat2=open(j,3);
    I=find(lon>=lon2,1); %find 1st entry
    J=find(lat>=lat2,1);
    if(I<=1 || J<=1) 
      disp('Failed to find an interval'); [j lon2 lat2 I J]
      error('Bomb out'); 
    end

    ratx=(lon(I)-lon2)/(lon(I)-lon(I-1));
    raty=(lat(J)-lat2)/(lat(J)-lat(J-1));
    %Check junks
    amp_max=max([amp0(I-1,J-1);amp0(I,J-1);amp0(I-1,J);amp0(I,J)]);
    amp_min=min([amp0(I-1,J-1);amp0(I,J-1);amp0(I-1,J);amp0(I,J)]);
    if(amp_max>100)
      if(amp_min>100 || amp_min<0)
        disp('All junks for amp:'); [i j I J]
        return;
      else %use min
        ap(i,1)=amp_min;
      end
    else
      tmp1=ratx*amp0(I-1,J-1)+(1-ratx)*amp0(I,J-1);
      tmp2=ratx*amp0(I-1,J)+(1-ratx)*amp0(I,J);
      ap(i,1)=tmp1*raty+tmp2*(1-raty);
    end %amp_max

    if(ap(i,1)>2)
      disp('Amp>2m');
      [i j amp_max amp_min]
    end

    %Make sure there is no jump in phase range
    wild(1:4)=[pha0(I-1,J-1); pha0(I,J-1); pha0(I,J); pha0(I-1,J)]; %counter-clockwise
    pha_max=max(wild);
    pha_min=min(wild);

    if(pha_max>370) %junk
      if(pha_min<0 ||  pha_min>370)
        disp('All junk phases:'); [i j I J]
        return;
      end
      ap(i,2)=pha_min;
    else
      for k=1:4
        if(abs(wild(k)-pha_max)>180)
          wild(k)=wild(k)+360;
          if(abs(wild(k)-pha_max)>180)
            disp('Phases still jump');
            [i j k I J]
            wild(1:4)
            return;
          end
        end
      end %for k

      tmp1=ratx*wild(1)+(1-ratx)*wild(2);
      tmp2=ratx*wild(4)+(1-ratx)*wild(3);
      ap(i,2)=tmp1*raty+tmp2*(1-raty);
    end %pha_max

    fprintf(fid,'%f %f\n',ap(i,1),ap(i,2));
  end %for j=1:npt
end %for i - freq's
fclose(fid);