Skip to content

Commit

Permalink
Merged BAFIM_flipchem from the flipchem branch
Browse files Browse the repository at this point in the history
  • Loading branch information
Ilkka Virtanen committed Jan 18, 2024
1 parent 72e854c commit 5407a27
Show file tree
Hide file tree
Showing 9 changed files with 1,618 additions and 33 deletions.
140 changes: 107 additions & 33 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,69 +1,143 @@
Bayesian Filtering Module (BAFIM)
# Bayesian Filtering Module (BAFIM)

[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4033903.svg)](https://doi.org/10.5281/zenodo.4033903)

BAFIM is an additional module to the GUISDAP incoherent scatter analysis tool ( https://www.eiscat.se/scientist/user-documentation/guisdap-9-0/ ). BAFIM replaces the standard IRI-based prior model with Bayesian filtering in time. Correlation priors are used in the prediction step of the filter to keep the predicted parameter profiles smooth in range direction in field-aligned measurements. The module can be used also with beam directions oblique to B and in analysis of remote data.
BAFIM is an additional module to the [GUISDAP](https://www.eiscat.se/scientist/user-documentation/guisdap-9-0/) incoherent scatter analysis tool. BAFIM replaces the standard IRI-based prior model with Bayesian filtering in time. Correlation priors are used in the prediction step of the filter to keep the predicted parameter profiles smooth in range direction in field-aligned measurements. The module can be used also with beam directions oblique to B and in analysis of remote data.

Reference: Virtanen, I.I, Tesfaw, H, Roininen, L., Lasanen, S., and Aikio, A., Bayesian filtering in incoherent scatter plasma parameter fits, manuscript submitted to JGR space physics, 2020.
The [Flipchem](https://github.com/amisr/flipchem) [![DOI](https://zenodo.org/badge/DOI/10.5281/10.5281/zenodo.3688853.svg)](https://doi.org/10.5281/zenodo.3688853) model can be used to support F region O+ ion fraction fits.

Ilkka Virtanen ([email protected]) 2020



#################################################################################################
INSTALLATION:
Ilkka Virtanen ([email protected]) 2024


## INSTALLATION

1. copy the files somewhere in your MATLAB path, for example the GUISDAP anal/ directory.

2. add the following lines in the end of the standard GUISDAP file anal/ionomodel.m:
2. add the following lines in the end of the standard GUISDAP function ionomodel in the file anal/ionomodel.m (before the last end statement)

apriorimodel=['apriorimodel_' iono_model];
if exist([apriorimodel '.m'])
[apriori,apriorierror] = feval(apriorimodel,apriori,apriorierror,heights,fit_altitude);
end



apriorimodel=['apriorimodel_' iono_model];
if exist([apriorimodel '.m'])
[apriori,apriorierror] = feval(apriorimodel,apriori,apriorierror,heights,fit_altitude);
end
## How to use the module?

There are two version of BAFIM available, the original BAFIM and the new BAFIM_flipchem. These are included as separate functions.


#################################################################################################
USAGE:
**Notice that user inputs are scaled differently in the present versions of BAFIM and BAFIM_flipchem!**

The scaling factors etc. mentioned below are explained in the reference paper.

### BAFIM (the original versoin)

The scaling factors etc. mentioned below are explained in the reference paper [Virtanen et al., 2021](https://doi.org/10.1029/2020JA028700).

To use BAFIM in GUISDAP fits:

1. Select 'bafim' as the ionospheric model in GUISDAP, i.e. write the following in the 'Special' box of the 'GUISDAP for dummies' window:

iono_model='bafim'
iono_model='bafim'

2. Set the lowest (hmin) and highest (hmax) altitude where each parameter is fitted, and scaling factors for process noise (st) and correlation length (sh) using the first four colums and six rows of the 'fit_altitude' array:

fit_altitude(1:6,1:4) = [ hmin_Ne , hmax_Ne , sh_Ne , st_Ne;
hmin_Ti , hmax_Ti , sh_Ti , st_Ti;
hmin_Tr , hmax_Tr , sh_Tr , st_Tr;
0 , 0 , 0 , 0 ;
hmin_Vi , hmax_Vi , sh_Vi , st_Vi;
hmin_Op , hmax_Op , sh_Op , st_Op ]

fit_altitude(1:6,1:4) = [ hmin_Ne , hmax_Ne , sh_Ne , st_Ne;
hmin_Ti , hmax_Ti , sh_Ti , st_Ti;
hmin_Tr , hmax_Tr , sh_Tr , st_Tr;
0 , 0 , 0 , 0 ;
hmin_Vi , hmax_Vi , sh_Vi , st_Vi;
hmin_Op , hmax_Op , sh_Op , st_Op ]
For example:

fit_altitude(1:6,1:4) = [ 0 , Inf , 0.1 , 2e11 ;
80 , Inf , 0.3 , 10 ;
103 , Inf , 0.3 , 0.05 ;
0 , 0 , 0 , 0 ;
80 , Inf , 0.2 , 2.5 ;
130 , 350 , 0.2 , 0.003 ]
fit_altitude(1:6,1:4) = [ 0 , Inf , 0.1 , 2e11 ;
80 , Inf , 0.3 , 10 ;
103 , Inf , 0.3 , 0.05 ;
0 , 0 , 0 , 0 ;
80 , Inf , 0.2 , 2.5 ;
130 , 350 , 0.2 , 0.003 ]

3. Disable the GUISDAP satellite check routine and set transmitter phase pushing to zero:
3. Tell the GUISDAP satellite check routine not to create gaps in the lag profiles and set transmitter phase pushing to zero:

a_phasepush=0
a_satch.do=0
a_phasepush=0
a_satch.cut=0

3. Start GUISDAP analysis as usual (hit 'GO')
4. Start GUISDAP analysis as usual (hit 'GO')




To run the RTS smoother, use the command:

bafim_smoother(<guisdap-output-dir>)
bafim_smoother(<guisdap-output-dir>)

After the smoothing step the Bayesian smoothing results are stored in variables r_param and r_error, as well as in r_param_smooth and r_error_smooth, while the Bayesian filtering outputs are in r_param_filter and r_error_filter. The results stored in the data files under r_param and r_error may be selected using the function bafim_select.




### BAFIM_flipchem

**The Flipchem package must be installed**. BAFIM_flipchem calls the Flipchem python module.

The scaling factors etc. mentioned below are explained in the reference paper [Virtanen et al., 2024](https://doi.org/10.22541/essoar.169945113.37163064/v1). **Notice that these are different from the original version**

To use BAFIM_flipchem in GUISDAP fits:

1. Select 'bafim' as the ionospheric model in GUISDAP, i.e. write the following in the 'Special' box of the 'GUISDAP for dummies' window:

iono_model='bafim_flipchem'

2. Set the lowest (hmin) and highest (hmax) altitude where each parameter is fitted, and scaling factors for process noise (st) and correlation length (sh) using the first four colums and six rows of the 'fit_altitude' array:

fit_altitude(1:6,1:4) = [ hmin_Ne , hmax_Ne , sh_Ne , st_Ne;
hmin_Ti , hmax_Ti , sh_Ti , st_Ti;
hmin_Tr , hmax_Tr , sh_Tr , st_Tr;
0 , 0 , 0 , 0 ;
hmin_Vi , hmax_Vi , sh_Vi , st_Vi;
hmin_Op , hmax_Op , sh_Op , st_Op ]
For example:

fit_altitude(1:6,1:4) = [ 0 , Inf , 0.05 , 2.5e11 ;
80 , Inf , 0.10 , 30 ;
97 , Inf , 0.10 , 0.03 ;
0 , 0 , 0 , 0 ;
80 , Inf , 0.05 , 2.5 ;
150 , 500 , 0.05 , 0.01 ]
**Notice that the numbers are different from those in the orignal BAFIM, because the length scales scale with time step duration and the chemistry model supports the composition fits.**

3. Tell the GUISDAP satellite check routine not to create gaps in the lag profiles and set transmitter phase pushing to zero:

a_phasepush=0
a_satch.cut=0

4. Start GUISDAP analysis as usual (hit 'GO')



To run the RTS smoother, use the command:

bafim_flipchem_smoother(<guisdap-output-dir>)


Four different results are stored in each data file after the smoothing step,

1. r_param_filter and r_error_filter contain the Bayesian filtering outputs, i.e. the parameters immediately after GUISDAP fit and Flipchem correction.
2. r_param_rcorr and r_error_rcorr contain the plasma parameters after applying the smoothing in range by correlation priors.
3. r_param_smooth and r_error_smooth are Bayesian smoothing results calculated from r_param_filter and r_error_filter. These results were shown in [Virtanen et al., 2021](https://doi.org/10.1029/2020JA028700).
4. r_param_rcorr_smooth and r_error_rcorr_smooth are Baeysian smoothing results calculated from r_param_rcorr and r_error_rcorr. These results are shown in [Virtanen et al., 2024](https://doi.org/10.22541/essoar.169945113.37163064/v1)

r_param_smooth and r_error smooth are written to r_param and r_error. However, this can be changed with the function bafim_flipchem_select to make any of these data easily readable by existing routines.

## References

Ilkka I. Virtanen, Habtamu W. Tesfaw, Anita T. Aikio, Roger Varney, Antti Kero, and Neethal Thomas, F1 region ion composition in Svalbard during the International Polar Year 2007-2008, Submitted to Journal of Geophysical Research: Space Physics, 2023, [https://doi.org/10.22541/essoar.169945113.37163064/v1](https://doi.org/10.22541/essoar.169945113.37163064/v1)

Ilkka I. Virtanen, Habtamu W. Tesfaw, Lassi Roininen, Sari Lasanen, and Anita Aikio, Bayesian filtering in incoherent scatter plasma parameter fits, Journal of Geophysical Research: Space Physics, 126 (3), 2021, [https://doi.org/10.1029/2020JA028700](https://doi.org/10.1029/2020JA028700)

Richards, P.G., Bilitza, D., Voglozin, D. Ion density calculator (IDC): A new efficient model of ionospheric ion densities. Radio Science, 45:1–11, 2010, [https://doi.org/10.1029/2009RS004332](https://doi.org/10.1029/2009RS004332).

Richards, P. G., Reexamination of ionospheric photochemistry. Journal of Geophysical Research: Space Physics, 116 (8), 1–15, [https://doi.org/10.1029/2011JA016613](https://doi.org/10.1029/2011JA016613)
43 changes: 43 additions & 0 deletions SS_flipchem.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
function SS = SS_flipchem(param,parinds,param_fit,covar,stdModErr,glat,glon,alt,fc)
%
% Sum-of-squares function to be minimized when adding flipchem information to the plasma parameter fits
%
% SS = SS_flipchem(param,parinds,param_fit,covar,stdModErr,date,glat,glon,alt)
%
% The flipchem python package must be installed (and the index files must be up to date)
%
% IV 2021
%


% first copy the original fit results
param_new = param_fit;

% then replace the new values given in param
param_new(parinds) = param;

% convert the parameters to real units
param_r = scaled_to_real(param_new);

% call the model and update the composition
ne = param_r(1);
te = param_r(2)*param_r(3);
ti = param_r(2);

% fc_comp = dirthe_flipchem(ne,te,ti,glat,glon,alt,fc);

outputs = fc.get_point(glat,glon,alt,ne,te,ti);

fc_comp = outputs{4}/ne;
% fc_comp = 1 - (outputs{5} + outputs{6})/ne;
% % take the scaled parameter vector, we will change only the composition
% param2 = param_new;
% param2(6) = outputsm{4}/param_r(1);

% % select the parameters that were actually fitted
% idd = diag(covar)~=0;

% SS = (param_fit(idd)-param_new(idd))*inv(covar(idd,idd))*(param_fit(idd)-param_new(idd))' + (modErr-priorModErr)^2/stdModErr^2;
SS = (param_fit(parinds)-param_new(parinds))*inv(covar(parinds,parinds))*(param_fit(parinds)-param_new(parinds))' + (param_new(6)-fc_comp)^2/stdModErr^2;

end
Loading

0 comments on commit 5407a27

Please sign in to comment.