%%%%%%%%%%%%%%%% TestModel time series %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
%%%%%%%%%%%%%%%%%%%%%%% Data Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MTimesteps=2000;
Lat=67; % degrees
Long=135; % degrees
MaxDeg=7; %% maximal degree of spherical harmonic expansion
modeltype='QC';
TimeStep = 10; %% years
cor_type = 'gaussian';
%%%%%%%%%%%%%%%%%%%%%%% OutPut File Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%
PrintYesNo=1;
formatstr='%+9.3f%+9.3f%+9.3f%+8.2f%+8.2f%9.3f%+8.2f%8.2f\n';
% FigYesNo=0;
% TimeScaleStr=sprintf('%dy',TimeStep*MTimesteps);

%%%%%%%%%%%%%%%%%%%%%%% Calculation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
latstep=15;
NumLat=4;


Lam2TimeStep = Memory2Lexic(MaxDeg,TimeStep); %% TimeStep to relative units
ZZ = StationaryGaussProcess(MTimesteps,Lam2TimeStep,cor_type); %% Base proc


    LatStr=sprintf('Lat%+.2f',Lat);
    DataFileName=sprintf('%s_Deg%d_%sCorr_Lat%+.2f_Long%+.2f_Tstep%dy',modeltype,MaxDeg,cor_type,Lat,Long,TimeStep);
    if PrintYesNo > 0        
        fid=fopen([DataFileName,'.txt'],'w');  
        fprintf(fid,'%7s%9s%9s%8s%8s%9s%10s%9s\n', 'X','Y','Z','I','D','F','VGPLat','VGPLong');    
    end
    NED_data=zeros(MTimesteps,3);
    IDF_data=zeros(MTimesteps,3);
    VGP_data=zeros(MTimesteps,2);

    Colat=90-Lat;
    [Harm,GradHarm] = HarmsAtPt(MaxDeg,Colat,Long);   %% Gradients of Harmonics
    
    [MeanM, SigmaM] = Model2Lexic(MaxDeg,modeltype);  %%  mean/sigma Parameters


    NLexicCff=size(Harm,1);
    for ii=1:MTimesteps
        %%make true process
        GaussCff=(ZZ(ii,:))'.*SigmaM + MeanM; %% gauss coefficients at time t_i
        for jj=1:NLexicCff
            NED_data(ii,:)=NED_data(ii,:)+GradHarm(jj,:)*GaussCff(jj);
        end

        IDF_data(ii,:)=NED2IDF(NED_data(ii,:));
        VGP_data(ii,:)=IncDec2VGP(IDF_data(ii,1),IDF_data(ii,2),Lat,Long);
        if PrintYesNo > 0
            fprintf(fid,formatstr, NED_data(ii,1),NED_data(ii,2),NED_data(ii,3), IDF_data(ii,1),IDF_data(ii,2),IDF_data(ii,3),VGP_data(ii,1),VGP_data(ii,2));
        end

    end
    
    if PrintYesNo > 0
         fclose(fid);
    end

clear

