clc
clear all
close all

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

%% init major tides
% -----------------
darwin = {'OM1','OM2','SA','SSA','MM','MF','MTM','MSQ','Q1','O1','P1','K1','2N2','N2','M2','S2','K2','M4'};

% 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

%% init all tides and admittance
% ------------------------------
% 256 tides
doodson ={'055.565','055.575','056.554','056.556','057.355','057.553','057.555','057.565','057.575',...
    '058.554','059.553','062.656','063.645','063.655','063.665','064.456','064.555','065.445','065.455',...
    '065.465','065.655','065.665','065.675','066.454','067.455','067.465','071.755','072.556','073.545',...
    '073.555','073.565','074.556','075.345','075.355','075.365','075.555','075.565','075.575','076.554',...
    '077.355','077.365','081.655','082.656','083.445','083.455','083.655','083.665','083.675','084.456',...
    '085.255','085.455','085.465','085.475','086.454','091.555','092.556','093.355','093.555','093.565',...
    '093.575','095.355','095.365','107.755','109.555','115.845','115.855','117.645','117.655','118.654',...
    '119.455','125.745','125.755','126.556','126.754','127.545','127.555','128.554','129.355','133.855',...
    '134.656','135.435','135.635','135.645','135.655','135.855','136.555','136.654','137.445','137.455',...
    '137.655','137.665','138.454','139.455','143.535','143.745','143.755','144.546','144.556','145.535',...
    '145.545','145.555','145.755','145.765','146.554','147.355','147.555','147.565','148.554','153.645',...
    '153.655','154.656','155.435','155.445','155.455','155.645','155.655','155.665','155.675','156.555',...
    '156.654','157.445','157.455','157.465','158.454','161.557','162.556','163.545','163.555','163.755',...
    '164.554','164.556','165.545','165.555','165.565','165.575','166.554','167.355','167.555','167.565',...
    '168.554','172.656','173.445','173.645','173.655','173.665','174.456','174.555','175.445','175.455',...
    '175.465','175.655','175.665','175.675','176.454','182.556','183.545','183.555','183.565','185.355',...
    '185.365','185.555','185.565','185.575','191.655','193.455','193.465','193.655','193.665','195.255',...
    '195.455','195.465','195.475','207.855','209.655','215.955','217.755','219.555','225.855','227.645',...
    '227.655','228.654','229.455','234.756','235.745','235.755','236.556','236.655','236.754','237.545',...
    '237.555','238.554','239.355','243.635','243.855','244.656','245.435','245.645','245.655','246.456',...
    '246.555','246.654','247.445','247.455','247.655','248.454','253.535','253.755','254.556','255.535',...
    '255.545','255.555','255.557','255.755','255.765','256.554','257.355','257.555','257.565','257.575',...
    '262.656','263.645','263.655','264.555','265.445','265.455','265.655','265.665','265.675','267.455',...
    '267.465','271.557','272.556','273.545','273.555','273.557','274.554','274.556','275.545','275.555',...
    '275.565','275.575','276.554','277.555','283.655','283.665','285.455','285.465','285.475','293.555',...
    '293.565','295.355','295.365','295.555','295.565','295.575','455.555'};

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

% interpolation matrix (256 tides -> 18 major tides)
admittance = load('data/admittance.linear.txt');


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

thetaf = doodsonMatrix * doodsonArguments(mjd); % 256 arguments for all tides
factorCos = admittance * cos(thetaf);  % interpolation coeff for 18 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}(1:maxDegree+1,1:maxDegree+1) + factorSin(i) * cnmSin{i}(1:maxDegree+1,1:maxDegree+1);
  snm = snm + factorCos(i) * snmCos{i}(1:maxDegree+1,1:maxDegree+1) + factorSin(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);
