%% init tide model
% ----------------
disp('init tide model');
path          = "../../models/EOT20_OCN/"; % path to .gfc files
basename      = "EOT20_OCN";        % filename prefix
doodsonMatrix = load (path+basename+'_002doodson.txt');
admittance    = load (path+basename+'_003admittance.txt');
fid           = fopen(path+basename+'_001fileList.txt');
fileList      = textscan(fid, '%s%s');
fclose(fid);

% read gfc files
for i=1:size(fileList{1}, 1)
  disp("read "+fileList{1}{i}+" and "+fileList{2}{i});
  anmCos{i} = readPotentialCoefficients(path+fileList{1}{i});
  anmSin{i} = readPotentialCoefficients(path+fileList{2}{i});
end

%% compute stokes coeffcients at time mjd
% ---------------------------------------
disp("compute stokes coeffcients at time mjd...");
mjd = 51544.5; % date for test computation

thetaf    = doodsonMatrix * doodsonArguments(mjd);
factorCos = admittance * cos(thetaf);
factorSin = admittance * sin(thetaf);

anm = zeros(size(anmCos{1}, 1), size(anmCos{1}, 2));
for i=1:size(fileList{1}, 1)
  anm = anm + factorCos(i) * anmCos{i} + factorSin(i) * anmSin{i};
end

% set degree one zero: Center of mass (CM)
anm(1:2, 1:2) = 0.0;

% anm contain now the stokes coeffcients of the ocean tides at time mjd
% Cnm coefficients of degree n and order m are stored in the lower triangle of anm,
% Snm coefficients of degree n and order m>0 are stored in the upper triangle of anm.

% ---------------------------------------

function dood=doodsonArguments(mjd)

t = (mjd-51544.5)/365250.0;

% compute gmst
era  = 2*pi * (mod(mjd-51544.5, 1) + 0.7790572732640 + 0.00273781191135448 * (mjd-51544.5));
gmst = era + (0.014506 + (4612.156534 + (1.3915817 + (-0.00000044 + (-0.000029956 + -0.0000000368*t)*t)*t)*t)*t) * pi/(180*3600);

dood = zeros(6,1);
dood(2,1) = (218.31664562999 + (4812678.81195750 + (-0.14663889 + ( 0.00185140 + -0.00015355*t)*t)*t)*t) * pi/180;
dood(3,1) = (280.46645016002 + ( 360007.69748806 + ( 0.03032222 + ( 0.00002000 + -0.00006532*t)*t)*t)*t) * pi/180;
dood(4,1) = ( 83.35324311998 + (  40690.13635250 + (-1.03217222 + (-0.01249168 +  0.00052655*t)*t)*t)*t) * pi/180;
dood(5,1) = (234.95544499000 + (  19341.36261972 + (-0.20756111 + (-0.00213942 +  0.00016501*t)*t)*t)*t) * pi/180;
dood(6,1) = (282.93734098001 + (     17.19457667 + ( 0.04568889 + (-0.00001776 + -0.00003323*t)*t)*t)*t) * pi/180;
dood(1,1) = gmst + pi - dood(2,1);
end

% ---------------------------------------

function anm=readPotentialCoefficients(filename)

fid = fopen(filename);

% read header
while 1
  line = fgetl(fid);
  if ~ischar(line), break, end
  if isempty(line), continue, end
  keyword = textscan(line,'%s',1);

  if(strcmp(keyword{1}, 'max_degree'))
    cells = textscan(line,'%s%d',1);
    maxDegree = cells{2};
    anm = zeros(maxDegree+1, maxDegree+1);
  end

  if(strcmp(keyword{1}, 'gfc'))
    cells = textscan(line,'%s%d%d%f%f');
    anm(cells{2}+1, cells{3}+1) = cells{4}; % cnm
    if cells{3}>0
      anm(cells{3}, cells{2}+1) = cells{5}; % snm
    end
  end
end

fclose(fid);
end

% ---------------------------------------
