%% This driver includes all the necessary definitions & finishes with call to main script
clear all

rad = pi/180;


% Following paths to directories with required sub-routines

warning off MATLAB:divideByZero

% Define vector datasets
options.filename_vector = char(...
    'dummy', ...
    'dummy',...
    'SWA_B_01.quiet', ...
    'dummy', ...
    'dummy');

% Define scalar datasets
options.filename_scalar = char(...
    'dummy', ...
    'dummy',...
    'SWA_F_01.quiet', ...
    'dummy',...
    'dummy');

% Define starting Euler values
params.Euler_start = [[0 0 0];
    [0 0 0];
    [11.8 -76.2 -12.6];
    [0 0 0];
    [0 0 0]];
params.Euler_weight = [0 0 0 0 0];
params.dt_Euler = [1e9 1e9 1e9 1e9 1e9];

% Dipole pole
params.igrf_dipole = [-29442.0,-1501.0, 4797.1];  

% Options related to GSM and SM co-ordinates
params.filename_GSM_SM_SHA = './GSM_SM_SHA_3D_IGRF12_2015.mat';
params.SM_i_e_flag = 'e+i';
params.GSM_i_e_flag = 'e+i';

options.filename_obs_data = 'dummy';

% read RC index
options.filename_Dst = '/nfs/wwwfiles/projects/oersted/models/RC/current/RC_1997-2022_augmented.dat';

% read Kp
options.filename_Kp = '/r2/solsys/indices/Kp_MJD_1998_QL.dat';

% file with E_m (merging electric field at magnetopause, for selection of polar data)
% options.filename_Em = 'D:\Matlab program\Geomagnetic field modeling\smallfit\tools\Em_index_2013_2019.mat';

% Data file for file to calculated QD co-ordinates
options.filename_apexsh = '/home/nio/m/tools/apexsh_1980-2025.txt';


% Define weighting options
options.robust = 'Huber';
%   options.robust = 'LS';
options_update.robust = 'Huber';
options.downweight_polar = 0;
params.w_obs_0 = 1.0;

% Define dataset options
options.N_Data_Frac = 5000;
options.N_min = 50;
options.both_stcam_CH = 0;

% Define data related parameters, variances and thresholds
params.sigma_csc_2 = 2.5^2;
params.sigma_ovh_2 = 2.5^2;
params.psi = 10/3600*rad;
params.chi = 10/3600*rad;

params.mu = 5;
params.PEJ_d_theta = 10;
params.PEJ_theta_0 = 14;
params.PC_variance = 0^2;

params.lat_gm_max = 55;
params.lat_gm_polar_cap = 90;
params.lat_gm_max_statistics = 55;
params.dDst_max = 2;
params.kp_max = 23;
params.Em_max = 10;

params.w_outlier = 0.1;
params.w_outlier = 0.01;


% Define factors needed for damping
params.a = 6371.2;
params.r_c = 3485;

% Choice of damping norms
options.damping_function = inline('(n+1).^2./(2.*n+1)', 'n');
options.damping_function_theta = inline('1-sin(theta).^3', 'theta');
options.spline_constrain_derivative = 3;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define specific model parameterization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options.N_Data_Step = 2;
options.add_data_polar_gap = 0;

model.maxdeg_MF = 15;

model.maxdeg_spline = 10;
model.maxdeg_ext = 2;


params.delta_t_q10 = 5e9;
params.delta_t_qs11 = 30e9;


% % Define default model choices
model.sp_order = 4;   % 6 in CHAOS

model.knots = [jd2000(2021, 1, 1) jd2000(2021, 12, 31)]';
model.knots = augknt(model.knots, model.sp_order);

model.epoch = 2021.25;

params.mu = 0;

%% DEFAULT PARAMETERS
% params.lambda_t_ll = .33;
% params.lambda_t_g10 = 10;
% params.lambda_t_1 = 10;
% params.lambda_t_2 = params.lambda_t_1;
%params.lambda_cr = 0.0;

% PV ==> Damping low latitude (Lambda 3 in CHAOS paper)
params.lambda_t_ll = 0.0;

% PV ==> Damping axial dipole coefficient (paragraph juste before section 3.2 in CHAOS paper)
params.lambda_t_g10 = 0;

% PV ==> Damping initial and end point (Lambda 2 in CHAOS paper)
params.lambda_t_1 = 0;
params.lambda_t_2 = params.lambda_t_1;

% PV ==> Damping crustal field (Lambda 0 in CHAOS paper)
params.lambda_cr = 0.0;




%%%%%%%%%%%%%%%%%%%%%%%%%
% Directory, Filenames  %
%%%%%%%%%%%%%%%%%%%%%%%%%

model.name = 'SWA';
model.header = [model.name, ', ', num2str(params.lambda_t_ll, '%e'), ', mu: ', num2str(params.mu, '%e'), ', w_obs_0: ', num2str(params.w_obs_0, '%e')];

options.filename_start_model = 'Swarm_11e_16.cof';


options.Max_iter = 6;

options.dir='./models/'; % Directory where wish to run from + where .cof file written to
options.dir_save='./models/'; % Directory where .mat files and other output is located

options.dir_update='./models/';
options.filename_update = [options.dir_update,'CHAOS_01j_14.mat'];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update : change options and params below to change details of model solved for in update %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
options.update = 0;

model_update.name=model.name;   % Take the model name as defined above for normal run
model_update.header=[model.name, ', '];

options_update.Max_iter=1;
params_update.lambda_t_ll = 10;
params_update.lambda_t_1  = 10;
params_update.lambda_t_2  = params_update.lambda_t_1 ;


diary off; diary 'Swarm_02_22.log';
model.name = 'Swarm_02a_22';

model.header = [model.name, ', ', num2str(params.lambda_t_ll, '%e'), ', mu: ', num2str(params.mu, '%e'), ', w_obs_0: ', num2str(params.w_obs_0, '%e')];
smallfit_noMPI_MSS;
