clear all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%
%            MONOPOL MODEL TIKHONOV REGULARIZED WITH HUBER WEIGHTS
%
%        -input: CHAOS4alpha vector and scalar data    
%        -input: MF7 vector/scalar synthetic data in polar regions
%                up til SH degree 60
%        -uses an icosahedron grid for monopole distribution
%
%        -output: radial component B_r
%        -Output: SH coefficients 
%
%
% Author:
% Magnus Danel Hammer
% Master Thesis project 2011, Copenhagen University
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                         CHAOS-4 vector data
disp('Load CHAOS-4 vector data')
tic
VEC      = load('/home/magdh/programCHAOS/data/CHAOS-4alpha-h-vector_NEC.dat');
%==========================================================================
%                         Data selected after time
%date = 3288.35342755;                       %1 januar 2009 -> 08:28:56
date    = 3288;
target2 = date;
wii     = find(VEC(:,1)>target2);
VEC     = VEC(wii,:);

t_vec           = VEC(:,1);                  %time [MD2000] in days (modified julian day)
r_vec           = VEC(:,2);                  %radius [km]
theta_vec       = VEC(:,3);                  %geographic co-latitude theta [deg]
phi_vec         = VEC(:,4);                  %geographic longitude phi [deg]
Bobs_r_vec      = VEC(:,5);                  %observed magnetic field, B_r [nT]                          
Bobs_theta_vec  = VEC(:,6);                  %observed magnetic field, B_theta [nT] 
Bobs_phi_vec    = VEC(:,7);                  %observed magnetic field, B_phi [nT] 
Bint_r_vec      = VEC(:,8);                  %internal part of model field, B_r [nT] 
Bint_theta_vec  = VEC(:,9);                  %internal part of model field, B_theta [nT] 
Bint_phi_vec    = VEC(:,10);                 %internal part of model field, B_phi [nT]
Bext_r_vec      = VEC(:,11);                 %external part of model field, B_r [nT]  
Bext_theta_vec  = VEC(:,12);                 %external part of model field, B_theta [nT] 
Bext_phi_vec    = VEC(:,13);                 %external part of model field, B_phi [nT]
%==========================================================================
%                   Data from VEC in vectors.
B_obs_vec = [Bobs_r_vec Bobs_theta_vec Bobs_phi_vec];                
B_ext_vec = [Bext_r_vec Bext_theta_vec Bext_phi_vec];                
B_int_vec = [Bint_r_vec Bint_theta_vec Bint_phi_vec];      %OBS! int contains both core and crustal fields!
%==========================================================================
%Create internal core field to degree n=15 from synth_value
%1)Get spline coefficients
load ('/home/magdh/programCHAOS/data/CHAOS-4alpha.mat')

%2)Construct new spline structure for SH model only up to degree N=15
N          = 15;
pp_N       = pp;
pp_N.dim   = N*(N+2);
coefs_tmp  = reshape(pp.coefs, [], pp.pieces, pp.order);
pp_N.coefs = reshape(coefs_tmp(1:N*(N+2),:,:), [], pp.order);
B_MF_vec   = synth_values(r_vec, theta_vec, phi_vec, pp_N, t_vec);

%Calculation of anomaly field by using B_MF_vec for core field
delB   = B_obs_vec(:,1) - (B_MF_vec(:,1) + B_ext_vec(:,1));   
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                         CHAOS-4 scalar data
disp('Load CHAOS-4 scalar data')
tic
SCA      = load('/home/magdh/programCHAOS/data/CHAOS-4alpha-h-scalar.dat');
%==========================================================================
%                           Data selected after time
%date    = 3288.68049931;                   %1 januar 2009 kl.16.19.55 
date    = 3288;
target2 = date;
vii     = find(SCA(:,1)>target2);
SCA     = SCA(vii,:);

t_sca     = SCA(:,1);                       %time [MD2000] in days (modified julian day)
r_sca     = SCA(:,2);                       %radius [km]
theta_sca = SCA(:,3);                       %geographic co-latitude theta [deg]
phi_sca   = SCA(:,4);                       %geographic longitude phi [deg]
Fobs_sca  = SCA(:,5);                       %observed magnetic field intensity F [nT]
Fint_sca  = SCA(:,6);                       %internal model prediction, F [nT]  %OBS! int contains both core and crustal fields!
Fext_sca  = SCA(:,7);                       %external model prediction, F [nT]
%==========================================================================
%Create internal core field to n=15 from synth_value
%1)Get spline coefficients
load ('/home/magdh/programCHAOS/data/CHAOS-4alpha.mat')

%2)Construct new spline structure for SH model only up to degree N=15
N          = 15;
pp_N       = pp;
pp_N.dim   = N*(N+2);
coefs_tmp  = reshape(pp.coefs, [], pp.pieces, pp.order);
pp_N.coefs = reshape(coefs_tmp(1:N*(N+2),:,:), [], pp.order);
B_MF_sca   = synth_values(r_sca, theta_sca, phi_sca, pp_N, t_sca);
F_MF_sca   = sqrt(B_MF_sca(:,1).^2 + B_MF_sca(:,2).^2 + B_MF_sca(:,3).^2);

%3)Calculate intensity anomaly
delF = Fobs_sca - (F_MF_sca + Fext_sca);
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                   ADD MF7 SCALAR DATA IN POLAR REGIONS
% 
%==========================================================================
disp('Generate synthetic MF7 polar scalar data')
tic
%North Pole: from 87 to 90 degrees
step    = 0.5;
theta_N = (0.5:step:6)';
phi_N   = (-180:step:180)';
[A1,B1] = ndgrid(theta_N,phi_N);
X1      = reshape(A1,1,[])';
Y1      = reshape(B1,1,[])';
r_N     = (6371.2+200).*ones(length(X1),1);

%South Pole: from -87 to -90 degrees
step    = 0.5;
theta_S = (174:step:180)';
phi_S   = (-180:step:180)';
[A2,B2] = ndgrid(theta_S,phi_S);
X2      = reshape(A2,1,[])';
Y2      = reshape(B2,1,[])';
r_S     = (6371.2+200).*ones(length(X2),1);
%==========================================================================
r_p     = [r_N; r_S];
theta_p = [X1; X2];
phi_p   = [Y1; Y2];
%==========================================================================
%==========================================================================
%load MF7 and CHAOS-4 main model
filename1 = '/home/magdh/programCHAOS/data/MF7.cof';
filename2 = '/home/magdh/programCHAOS/data/CHAOS-2s_alpha_04.cof'
%filename3 = '/home/magdh/programCHAOS/data/CHAOS-4alpha.mat';
%==========================================================================
%create synthetic scalar values from MF7 in CHAOS-4 scalarlocations - polar synthetic data
[m_i_sca, epoch, maxdeg_MF, maxdeg_SV, r_ref, m_ext, maxdeg_SV2, maxdeg_SV3, model_header] = read_model_3(filename1);

%maximum SH degree and order N=60; 
deg_pol = 60;
m_i_sca = m_i_sca(1:deg_pol*(deg_pol+2));
m_MF_p  = m_i_sca;
delB_p  = synth_values(r_p, theta_p, phi_p, m_MF_p);   
%==========================================================================
%linearization by CHAOS-2 main model - polar synth data
[m_MF_2, epoch, maxdeg_MF, maxdeg_SV, r_ref, m_ext, maxdeg_SV2, maxdeg_SV3, model_header] = read_model_3(filename2);

B_MF_p = synth_values(r_p, theta_p, theta_p, m_MF_2(1:255)); %N=15 (15*17)
F_MF_p = sqrt(B_MF_p(:,1).^2 + B_MF_p(:,2).^2 + B_MF_p(:,3).^2);
delF_p = (delB_p(:,1).*B_MF_p(:,1) + delB_p(:,2).*B_MF_p(:,2) + delB_p(:,3).*B_MF_p(:,3))./F_MF_p;
%==========================================================================
%linearization by CHAOS-4 main model - polar synth data
%load ('/home/magdh/programCHAOS/data/CHAOS-4alpha.mat')
% t_p        = [1:1:length(r_p)]';
% N          = 15;
% pp_N       = pp;
% pp_N.dim   = N*(N+2);
% coefs_tmp  = reshape(pp.coefs, [], pp.pieces, pp.order);
% pp_N.coefs = reshape(coefs_tmp(1:N*(N+2),:,:), [], pp.order);
% B_MF_p     = synth_values(r_p, theta_p, phi_p, pp_N, t_p);
% F_MF_p     = sqrt(B_MF_p(:,1).^2 + B_MF_p(:,2).^2 + B_MF_p(:,3).^2);
 
% %scalar anomaly field linearization - polar synth data
% delF_p = (delB_p(:,1).*B_MF_p(:,1) + delB_p(:,2).*B_MF_p(:,2) + delB_p(:,3).*B_MF_p(:,3))./F_MF_p;
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc

%totale antal data punkter
r_data = [r_vec; r_sca];
r_tot  = [r_vec; r_sca; r_p];


tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                    DISTRIBUTION OF POINTS SOURCES
%
%==========================================================================
%Create an 'iso' grid on the sphere using sphere_tri
disp('Distribution of point sources')
maxlevel = 5;
a        = 6371.2; 
winding  = 1;
FV       = sphere_tri('ico',maxlevel,a,winding)
%==========================================================================
%calculate midpoints of all triangles. 
MP    = zeros(size(FV.faces,1), 3);
for n = 1:size(FV.faces,1)
     MP(n,:) = [sum(FV.vertices(FV.faces(n,:),1)) sum(FV.vertices(FV.faces(n,:),2)) sum(FV.vertices(FV.faces(n,:),3))]/3;
end
%==========================================================================
%Using centerpoints and/or vertices in grid
Y=[FV.vertices;MP];
%Y=MP;
%Y=FV.vertices;
MONO=Y;
%==========================================================================
%                           COORDINATES
%==========================================================================
%Sourcepunkterne is given by (s,theta_s,phi_s) from MONO. 
%Coordinates for satellite data is given by (r,theta,phi).
%==========================================================================
%monopoldistance in km (distance between monopoles):
a       = 6371.2;                           %Earth radius
mon_km  = sqrt(4*pi*a^2/(length(Y)))

%monopoldistance in degrees (1degree=111.3km):
mon_g   = sqrt(4*pi*a^2/(length(Y)))/111.3
%==========================================================================
rad     = pi/180;                           %omregning til radianer
deg     = 180/pi;                           %omregning til grader
%==========================================================================
depth   = 200;                              %depth of monopoles  
s       = (6371.2-depth)*ones(length(Y),1);  
theta_s = acosd(Y(:,3)./sqrt(Y(:,1).^2+Y(:,2).^2+Y(:,3).^2));
phi_s   = atan2(Y(:,1),Y(:,2))*deg; 
%========================================================================== 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc

disp('save data and model locations and data observations')
tic
i=1;
save(sprintf('/home/magdh/NewCalc/CHAOS_dataT_%1.1d.mat',i), '-v7.3','B_MF_sca','F_MF_sca','B_MF_p','F_MF_p','r_vec','r_sca','r_p','theta_vec','theta_sca','theta_p','phi_vec','phi_sca','phi_p','delB','delF','delF_p','s','theta_s','phi_s','delF_p','r_p','theta_p','phi_p');
toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%clear variables
clear('vii','wii','pp','pp_N','N','g','date','coefs_tmp','dg','target2')
clear('B_obs_vec','B_ext_vec','B_int_vec','Bext_r_vec','Bext_theta_vec','Bext_phi_vec','Bint_r_vec','Bint_theta_vec','Bint_phi_vec','VEC')
clear('Bobs_r_vec','Bobs_theta_vec','Bobs_phi_vec','B_MF_vec')
clear('Fobs_sca','Fint_sca','Fext_sca','SCA')
clear('MP','maxlevel','winding','FV')
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                               IRLS BEGIN
%
%==========================================================================
%generates initial modelvector and weighting matrix for vector and scalar 
%data
disp('Generate initial modelvector and weighting matrix')
tic
m_r = zeros(length(phi_s),1);
m   = zeros(length(phi_s),1);
w_r = ones(length(phi_vec),1);
w_F = ones(length(phi_sca),1);
w_P = ones(length(phi_p),1);
toc
%==========================================================================
%number of iterations                           
Max_iter = 5;                   
for iter = 1:Max_iter
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                        VECTOR GREEN MATRIX 
%                          - r dependence
%==========================================================================
disp('Calculation of vector Green matrix, iter')
tic
delB_r_synth = zeros(length(phi_vec),1);
GtWG    = zeros(length(phi_s),length(phi_s));
GtWy    = zeros(length(phi_s), 1);
sin_t   = sind(theta_vec);
cos_t   = cosd(theta_vec);
sin_t_s = sind(theta_s);
cos_t_s = cosd(theta_s);
N_step  = 1000;                    
hwait=waitbar(0,'Please wait - vector Green matrix');
for k = 1:N_step:length(phi_vec);
    waitbar(i/length(phi_vec),hwait);
    G_r     = zeros(N_step,length(phi_s));
    n = 0;
    k_max = min(k+N_step, length(phi_vec));
    for i = k:k_max
        d_vec    = sqrt(s.^2+r_vec(i).^2-(2.*r_vec(i).*s.*(sin_t(i).*sin_t_s.*cosd(phi_vec(i)-phi_s)+ cos_t(i).*cos_t_s)));
        n        = n+1;
        G_r(n,:) = ((r_vec(i) - s.*(sin_t(i).*sin_t_s.*cosd(phi_vec(i)-phi_s)+ cos_t(i).*cos_t_s)  )./ (d_vec.^3))';
    end
    GtWG = GtWG + G_r(1:n,:)'*(repmat(w_r(k:k_max),1,length(phi_s)).*G_r(1:n,:));
    GtWy = GtWy + G_r(1:n,:)'*(w_r(k:k_max).*delB(k:k_max,1));
    delB_r_synth(k:k_max) = G_r(1:n,:)*m;
end
close(hwait)
toc
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                         SCALAR GREEN MATRIX - DATA
%
%==========================================================================
disp('Calculation of scalar Green matrix - data')
tic
delF_synth = zeros(length(phi_sca),1);
AtWA    = zeros(length(phi_s),length(phi_s));
AtWy    = zeros(length(phi_s), 1);
sin_ts  = sind(theta_sca);
cos_ts  = cosd(theta_sca);
sin_t_s = sind(theta_s);
cos_t_s = cosd(theta_s);
N_step  = 1000; 
hwait=waitbar(0,'Please wait - scalar Green matrix (data)');
for k = 1:N_step:length(phi_sca);
    waitbar(i/length(phi_sca),hwait);
    A_r  = zeros(N_step,length(phi_s));
    A_t  = zeros(N_step,length(phi_s));
    A_p  = zeros(N_step,length(phi_s));
    A_F  = zeros(N_step,length(phi_s));
    n = 0;
    k_max = min(k+N_step, length(phi_sca));
    for i = k:k_max
            n       = n+1;
        d_sca       = sqrt(s.^2+r_sca(i).^2-(2.*r_sca(i).*s.*(sin_ts(i).*sin_t_s.*cosd(phi_sca(i)-phi_s)+ cos_ts(i).*cos_t_s))).^3;    
        A_r(n,:) = (( r_sca(i) - s.*(sin_ts(i).*sin_t_s.*cosd(phi_sca(i)-phi_s)+ cos_ts(i).*cos_t_s)  )./(d_sca))';
        A_t(n,:) = (( s.*(sin_ts(i).*cos_t_s - cos_ts(i).*sin_t_s.*cosd(phi_sca(i)-phi_s))            )./(d_sca))';
        A_p(n,:) = (( s.*sin_t_s.*sind(phi_sca(i)-phi_s)                                              )./(d_sca))';       
        A_F(n,:) = ( A_r(n,:).*B_MF_sca(i,1) + A_t(n,:).*B_MF_sca(i,2) +A_p(n,:).*B_MF_sca(i,3)       )./F_MF_sca(i);
    end 
     AtWA                = AtWA + A_F(1:n,:)'*(repmat(w_F(k:k_max),1,length(phi_s)).*A_F(1:n,:));
     AtWy                = AtWy + A_F(1:n,:)'*(w_F(k:k_max).*delF(k:k_max));
     delF_synth(k:k_max) = A_F(1:n,:)*m;    
end
close(hwait)
toc
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                    SCALAR GREEN MATRIX -POLAR SYNTH DATA
%
%==========================================================================
disp('Calculation of scalar Green matrix - synth')
tic
delF_synth_p = zeros(length(phi_p),1);
AtWA_p    = zeros(length(phi_s),length(phi_s));
AtWy_p    = zeros(length(phi_s), 1);
sin_ts_p  = sind(theta_p);
cos_ts_p  = cosd(theta_p);
sin_t_s   = sind(theta_s);
cos_t_s   = cosd(theta_s);
N_step  = 1000; 
hwait=waitbar(0,'Please wait - scalar Green matrix (synth)');
for k = 1:N_step:length(phi_p);
    waitbar(i/length(phi_p),hwait);
    A_r_p  = zeros(N_step,length(phi_s));
    A_t_p  = zeros(N_step,length(phi_s));
    A_p_p  = zeros(N_step,length(phi_s));
    A_F_p  = zeros(N_step,length(phi_s));
    n = 0;
    k_max = min(k+N_step, length(phi_p));
    for i = k:k_max
            n       = n+1;
         d_p        = sqrt(s.^2+r_p(i).^2-(2.*r_p(i).*s.*(sin_ts_p(i).*sin_t_s.*cosd(phi_p(i)-phi_s)+ cos_ts_p(i).*cos_t_s))).^3;   
        A_r_p(n,:) = (( r_p(i) - s.*(sin_ts_p(i).*sin_t_s.*cosd(phi_p(i)-phi_s)+ cos_ts_p(i).*cos_t_s)  )./(d_p))';
        A_t_p(n,:) = (( s.*(sin_ts_p(i).*cos_t_s - cos_ts_p(i).*sin_t_s.*cosd(phi_p(i)-phi_s))          )./(d_p))';
        A_p_p(n,:) = (( s.*sin_t_s.*sind(phi_p(i)-phi_s)                                                )./(d_p))';       
        A_F_p(n,:) = ( A_r_p(n,:).*B_MF_p(i,1) + A_t_p(n,:).*B_MF_p(i,2) +A_p_p(n,:).*B_MF_p(i,3)       )./F_MF_p(i);
    end 
    AtWA_p                = AtWA_p + A_F_p(1:n,:)'*(repmat(w_P(k:k_max),1,length(phi_s)).*A_F_p(1:n,:));
    AtWy_p                = AtWy_p + A_F_p(1:n,:)'*(w_P(k:k_max).*delF_p(k:k_max));
    delF_synth_p(k:k_max) = A_F_p(1:n,:)*m;    
end
close(hwait)
toc
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                        COMBINING DESIGN MATRICES
% 
%==========================================================================
%combine design matrix
disp('Combination of design matrix')
tic
GtWG = GtWG + AtWA + AtWA_p;
GtWy = GtWy + AtWy + AtWy_p;   
toc

%find unconstrained solution (i.e. mean(m_r)=0 not fulfilled) with regularization 
disp('Solution by regularisation')
lambda = 1e-13;
%lambda = 0;
tic
m_r = (GtWG + lambda.*eye(size(GtWG,1)))\GtWy; 
toc

%constrained solution mean(m)=0
disp('Constrained solution')
tic
L = ones(length(m_r),1);
m = m_r - (GtWG)^(-1)*L*(L'*m_r)*(L'*(GtWG)^(-1)*L)^(-1); 
toc
%tjek solution
mean(m)
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                           HUBER WEIGHTING
% 
%==========================================================================
%VECTOR WEIGHTING
disp('Weighting of vector data')
tic
%residuals between observations and model predictions
res_r = (delB(:,1) - delB_r_synth);

%weighted mean of residuals
MEAN_r = sum(w_r.*res_r)/sum(w_r)

%root mean square
rms_r = sqrt(sum(w_r.*(delB(:,1) - delB_r_synth).^2)/sum(w_r))

%standard deviation
sigma_r = rms_r;

%Huber weights for w_r
c = 1.5;
for l=1:length(delB(:,1))
eps_r(l) = (abs(res_r(l))/sigma_r);
k(l)     = c/eps_r(l);
    if k(l) > abs(eps_r(l))
        w_r(l) = 1;
    else 
        w_r(l) = c/eps_r(l); 
    end
end
toc
%==========================================================================
%SCALAR WEIGHTING
disp('Weighting of scalar data - data')
tic
%residuals between observations and model predictions
res_F = (delF - delF_synth);

%weighted mean of residuals
MEAN_F = sum(w_F.*res_F)/sum(w_F)

%root mean square
rms_F = sqrt(sum(w_F.*(delF - delF_synth).^2)/sum(w_F))

%standard deviation
sigma_F = rms_F;

%Huber weights for W_F
c = 1.5;
for t=1:length(delF)
eps_F(t) = (abs(res_F(t))/sigma_F);
k(t)     = c/eps_F(t);
    if k(t) > abs(eps_F(t))
        w_F(t) = 1;
    else 
        w_F(t) = c/eps_F(t); 
    end
end
toc
%==========================================================================
%SCALAR WEIGHTING
disp('Weighting of scalar data - synth')
tic
%residuals between observations and model predictions
res_F_p = (delF_p - delF_synth_p);

%weighted mean of residuals
MEAN_F_p = sum(w_P.*res_F_p)/sum(w_P)

%root mean square
rms_F_p = sqrt(sum(w_P.*(delF_p - delF_synth_p).^2)/sum(w_P))

%standard deviation
sigma_F_p = rms_F_p;

%Huber weights for W_F
c = 1.5;
for t=1:length(delF_p)
eps_F_p(t) = (abs(res_F_p(t))/sigma_F_p);
k_p(t)     = c/eps_F_p(t);
    if k_p(t) > abs(eps_F_p(t))
        w_P(t) = 1;
    else 
        w_P(t) = c/eps_F_p(t); 
    end
end
toc
%==========================================================================
%save after each iteration
disp('save after each iteration')
tic
%save(sprintf('GreenMatrices_%1.1d.mat', iter), '-v7.3','GtWG','GtWy','AtWA','AtWy','AtWA_p','AtWy_p');
save(sprintf('/home/magdh/NewCalc/ModelValues30722_13_%1.1d.mat', iter), '-v7.3','w_F','w_r','w_P','m_r','res_r','MEAN_r','rms_r','MEAN_F','res_F','rms_F','res_F_p','rms_F_p','MEAN_F_p','m');
toc
end
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                               IRLS END
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%effektive brug af målinger efter vægtning
sum(w_r)/length(w_r)
sum(w_F)/length(w_F)
sum(w_P)/length(w_P)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%==========================================================================
%                   MONOPOLE MODEL SH GAUSS COEFFICIENTS
% 
%   -uses program: GaussCoef.m (Magnus Danel Hammer)
%   -uses program: powerspec.m (Nils Olsen)
%==========================================================================
N=180;                                                  %maximum SH degree and order of model
[g,power] = GaussCoef(s,theta_s,phi_s,N,m);             %calculates coefficients and power  
save(sprintf('/home/magdh/NewCalc/GaussCoef_T30722_13_%1.1d.mat', iter), '-v7.3', 'g','power');
%==========================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

