%This is sample main script to help users understanding our model more easily

% POS is a array which contains all sites 'positions,
% POS first column stores all sites' longitudes,
% and second column must be all sites' latitudes.
clear all;
clc;
POS = [  -178.23 ,  35.85;     % first site: lon -178.23, and lat  35.85
          152.32 , -22.56];    % second site: lon 152.32, and lat -22.56

% Generate Ts-Tm function coefficients at given positions
%  Ts-Tm function: 
%   f(1)*X(:,1)+f(2)+f(3)*cos(X(:,2)*2*pi/365.25)...
%  +f(4)*sin(X(:,2)*2*pi/365.25)+f(5)*cos(X(:,2)*4*pi/365.25)...
%  +f(6)*sin(X(:,2)*4*pi/365.25)+f(7)*cos(X(:,3)*2*pi/24)...
%  +f(8)*sin(X(:,3)*2*pi/24),
%   where f(1) slope constant,
%         f(2) intercpet constant,
%         f(3) and f(4) for annual variations terms,
%         f(5) and f(6) for semiannual variations terms,
%         f(7) and f(8) for diurnal variations terms,
%   and   X is the function input array,
%         X(:,1), first column in array X, surface temperature input series
%         X(:,2), UTC day of year(doy),
%         X(:,3), UTC hour of day .
Coeffs = Generate_Ts_Tm_coeffs(POS);


%%%%% Test with simulation data 
% Surface temperature, doy and hour inputs
INPUT = [ 280.5 , 100,  23 ;     %surface temperature 280.5 K, at UTC 23:00, doy 100
       295.2,  365,  10.5];   %surface temperature 295.2 K, at UTC 10:30, doy 365

% Ts-Tm function coefficients at the first site
C = Coeffs(1,:);

% Calculate Tm at the selected site
modelfun=@(f,X)f(1)*X(:,1)+f(2)...
    +f(3)*cos(X(:,2)*2*pi/365.25)+f(4)*sin(X(:,2)*2*pi/365.25)...
    +f(5)*cos(X(:,2)*4*pi/365.25)+f(6)*sin(X(:,2)*4*pi/365.25)...
    +f(7)*cos(X(:,3)*2*pi/24)+f(8)*sin(X(:,3)*2*pi/24);
TM = modelfun(C,INPUT);


