clc
clear all
close all

%% Constants
maxDegree = 180;
mjd = 51544.5; % date for test computation

%% init major tides
% -----------------
darwin = {'om1', 'om2', 'sa', 'ssa', 'mm', 'mf', 'mtm', 'msq', ...
          'q1', 'o1', 'p1', 's1', 'k1', 'j1', ...
          'eps2', '2n2', 'mu2', 'n2', 'nu2', 'm2', 'la2', 'l2', 't2', 's2', 'r2', 'k2', ...
          'm3', 'n4', 'mn4', 'm4', 'ms4', 's4', 'm6', 'm8'};

% read stokes coefficients
fileName = '../fes2014b/fes2014b_n180_version20170520';
for i=1:size(darwin,2)
  disp(strcat('read:', fileName,'.',darwin{i},'.gfc'));
  [cnmCos{i} snmCos{i} GM R] = readPotentialCoefficients( strcat(fileName,'.',darwin{i},'.cos.gfc') );
  [cnmSin{i} snmSin{i} GM R] = readPotentialCoefficients( strcat(fileName,'.',darwin{i},'.sin.gfc') );
end

%% init all tides and admittance
% ------------------------------
doodson = textscan(fopen('../admittance/admittanceTides.txt'), '%q', 'CommentStyle', '#');

% (361 x 6) matrix with doodson multipliers converted from the doodson strings
doodsonMatrix = zeros(size(doodson{1},1), 6);
for i=1:size(doodson{1},1)
  doodsonMatrix(i,1) = hex2dec(doodson{1}{i}(1))-0;
  doodsonMatrix(i,2) = hex2dec(doodson{1}{i}(2))-5;
  doodsonMatrix(i,3) = hex2dec(doodson{1}{i}(3))-5;
  doodsonMatrix(i,4) = hex2dec(doodson{1}{i}(5))-5;
  doodsonMatrix(i,5) = hex2dec(doodson{1}{i}(6))-5;
  doodsonMatrix(i,6) = hex2dec(doodson{1}{i}(7))-5;
end

% admittance interpolation matrix (361 tides -> 34 major tides)
admittance = load('../admittance/fes2014b_admittance_linear_linear.txt');


%% compute stokes coeffcients at time mjd
% ---------------------------------------
disp('compute stokes coeffcients at time mjd...');

thetaf = doodsonMatrix * doodsonArguments(mjd); % arguments for all tides
factorCos = admittance * cos(thetaf);           % interpolation coeff for the major tides
factorSin = admittance * sin(thetaf);

% sum up all major tides
cnm = zeros(maxDegree+1, maxDegree+1);
snm = zeros(maxDegree+1, maxDegree+1);
for i=1:size(darwin,2)
  cnm = cnm + factorCos(i) * cnmCos{i} + factorSin(i) * cnmSin{i};
  snm = snm + factorCos(i) * snmCos{i} + factorSin(i) * snmSin{i};
end
% cnm,snm contain now the stokes coeffcients of the ocean tides at time mjd


%% compute water height at 1x1 grid on sphere with radius R
% ---------------------------------------------------------
disp('compute water height at 1x1 grid on sphere with radius R...');
lambda = pi/180 * [-179.5:1:179.5];
phi    = pi/180 * [-89.5:1:89.5];
plotWaterHeight(lambda, phi, cnm, snm, GM, R);
