clc
clear all
close all

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

%% init major tides
% -----------------
darwin{1}  = 'OM1';  doodson{1}  = '055.565';
darwin{2}  = 'OM2';  doodson{2}  = '055.575';
darwin{3}  = 'SA';   doodson{3}  = '056.554';
darwin{4}  = 'SSA';  doodson{4}  = '057.555';
darwin{5}  = 'MM';   doodson{5}  = '065.455';
darwin{6}  = 'MF';   doodson{6}  = '075.555';
darwin{7}  = 'MTM';  doodson{7}  = '085.455';
darwin{8}  = 'MSQ';  doodson{8}  = '093.555';
darwin{9}  = 'Q1';   doodson{9}  = '135.655';
darwin{10} = 'O1';   doodson{10} = '145.555';
darwin{11} = 'P1';   doodson{11} = '163.555';
darwin{12} = 'K1';   doodson{12} = '165.555';
darwin{13} = '2N2';  doodson{13} = '235.755';
darwin{14} = 'N2';   doodson{14} = '245.655';
darwin{15} = 'M2';   doodson{15} = '255.555';
darwin{16} = 'S2';   doodson{16} = '273.555';
darwin{17} = 'K2';   doodson{17} = '275.555';
darwin{18} = 'M4';   doodson{18} = '455.555';

% read stokes coeffcients
fileName = 'eot11a/eot11a';
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

% 18x6 matrix with doodson multipliers
doodsonMatrix = zeros(size(doodson,2), 6);
for i=1:size(doodson,2)
  doodsonMatrix(i,1) = doodson{i}(1)-'0';
  doodsonMatrix(i,2) = doodson{i}(2)-'5';
  doodsonMatrix(i,3) = doodson{i}(3)-'5';
  doodsonMatrix(i,4) = doodson{i}(5)-'5';
  doodsonMatrix(i,5) = doodson{i}(6)-'5';
  doodsonMatrix(i,6) = doodson{i}(7)-'5';
end


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

thetaf = doodsonMatrix * doodsonArguments(mjd);

cnm = zeros(maxDegree+1, maxDegree+1);
snm = zeros(maxDegree+1, maxDegree+1);
for i=1:size(doodsonMatrix,1)
  cnm = cnm + cos(thetaf(i)) * cnmCos{i}(1:maxDegree+1,1:maxDegree+1) + sin(thetaf(i)) * cnmSin{i}(1:maxDegree+1,1:maxDegree+1);
  snm = snm + cos(thetaf(i)) * snmCos{i}(1:maxDegree+1,1:maxDegree+1) + sin(thetaf(i)) * snmSin{i}(1:maxDegree+1,1:maxDegree+1);
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);
