国际参考电离层2020模型Matlab（IRI-2020）
国际参考电离层 （IRI） 是一个国际项目 由空间研究委员会 （COSPAR） 和国际 无线电科学联盟 （URSI）。这些组织组成了一个工作组 （在 60 年代后期产生了电离层的经验标准模型， 基于所有可用数据源 （几个稳步改进的版本 模型已经发布。电子密度、电子温度、离子温度、 离子组成（O+、H+、He+、N+、NO+、O+2、簇离子）、赤道垂直离子漂移、垂直电离层电子含量（vTEC;
·
一、IRI模型
国际参考电离层 （IRI） 是一个国际项目 由空间研究委员会 （COSPAR）和国际无线电科学联盟 （URSI）。该组织组成了一个工作组 （ 60 年代后期至今）产生了电离层的经验标准模型，最新版本是IRI-2020。
对于给定的地点、时间和日期，IRI 模型基本参数如下：
输出参数：
电子密度、电子温度、离子温度、 离子组成（O+、H+、He+、N+、NO+、O+2、簇离子）、赤道垂直离子漂移、垂直电离层电子含量（vTEC;用户可以选择沿电子密度分布进行积分的起点和终点高度）、F1 层的出现概率 spread-F 和偶发E 的极光边界，电离层风暴对 F 和 E 峰值密度的影响，以及等离子体频率与电子回旋频率的比率。
电子密度、电子温度、离子温度、 离子组成（O+、H+、He+、N+、NO+、O+2、簇离子）、赤道垂直离子漂移、垂直电离层电子含量（vTEC;用户可以选择沿电子密度分布进行积分的起点和终点高度）、F1 层的出现概率 spread-F 和偶发E 的极光边界，电离层风暴对 F 和 E 峰值密度的影响，以及等离子体频率与电子回旋频率的比率。
输入参数：
必需：
太阳指数（F10.7 每日、81 天和 12 个月运行平均值;太阳黑子数 12 个月运行平均值）
电离层指数（基于电离探测仪的 IG 指数 12 个月运行平均值）
磁指数（3 小时 ap、每日 ap）
可选输入参数：
用户可以提供许多输入参数，IRI 配置文件将根据以下输入参数进行调整：
F2 峰高 （hmF2/km） 或传播因子 M（3000）F2 与 hmF2
F2 峰等离子体频率 （foF2/MHz） 或电子密度 （NmF2/m-3）
底部轮廓参数 B0/km（厚度）和 B1（形状）
F1 等离子体频率 （foF1/MHz） 或电子密度 （NmF1/m-3）
E 峰高 （hmE/km）
E 峰等离子体频率 （foE/MHz） 或电子密度 （NmE/m-3）
必需：
太阳指数（F10.7 每日、81 天和 12 个月运行平均值;太阳黑子数 12 个月运行平均值）
电离层指数（基于电离探测仪的 IG 指数 12 个月运行平均值）
磁指数（3 小时 ap、每日 ap）
可选输入参数：
用户可以提供许多输入参数，IRI 配置文件将根据以下输入参数进行调整：
F2 峰高 （hmF2/km） 或传播因子 M（3000）F2 与 hmF2
F2 峰等离子体频率 （foF2/MHz） 或电子密度 （NmF2/m-3）
底部轮廓参数 B0/km（厚度）和 B1（形状）
F1 等离子体频率 （foF1/MHz） 或电子密度 （NmF1/m-3）
E 峰高 （hmE/km）
E 峰等离子体频率 （foE/MHz） 或电子密度 （NmE/m-3）
高度范围：
电子密度：白天：65 - 30,000 公里，夜间：80 - 30,000 公里
电子和离子温度：60 - 2,000 公里（IRI-95 选项：60 - 3,000 公里）
离子组成：75-2,000 公里（DS95/DY85 选项：80-2,000 公里）
电子密度：白天：65 - 30,000 公里，夜间：80 - 30,000 公里
电子和离子温度：60 - 2,000 公里（IRI-95 选项：60 - 3,000 公里）
离子组成：75-2,000 公里（DS95/DY85 选项：80-2,000 公里）
现状：
官方提供了使用Fortran代码，Fortran计算效率高但是可视化效果不佳，使用MATLAB用于科学分析更加简便，但目前IRI官方网站提供的MATLAB版本是实时读取IRI网页数据并进行处理保存，因此需要联网使用，单机无法运行，不适合二次开发。最佳解决方案是mex把IRI模型的Fortran源代码封装成为matlab可以调用的函数。
二、接口函数
mex_iri.F
#include "fintrf.h"
C Gateway routine
subroutine mexFunction(nlhs, plhs, nrhs, prhs)
c MATLAB mex wrapper for the IRI 2020 model
c
c Accepts as input a size [9,n] array containing inputs for
c the NeQuick model. Each column should have the following:
c
C h: height (km)
C alat: gg. latitude (degrees N)
C along: gg. longitude (degrees E)
C mth: month (1 .. 12)
C flx: 10.7 cm solar radio flux (flux units)
C UT: Universal Time (hours)
c
c Must be run from the folder containing the ccirxx.asc
c and modip.asc files
c
C Output is size (n,1) electron density in units of m^-3
C (electrons per cubic meter)
c
C Declarations
C mexFunction arguments:
mwPointer plhs(*), prhs(*)
integer nlhs, nrhs
C Function declarations:
mwPointer mxGetPr
mwPointer mxCreateDoubleMatrix
integer mxIsNumeric
mwPointer mxGetM, mxGetN
integer*4 mexPrintf
C Local Variables
mwPointer in_ptr, outf_ptr, oar_ptr, uin_ptr
real*8 ijf(50), input(9), uparam(16)
real*4 outf(20,1000),oar(100)
LOGICAL*4 jf(50)
integer nummax, numhei, tnmax
real*4 heibeg, heiend, heistp
logical UINPUT
C Array information:
mwPointer mrows, ncols, umrows, uncols
mwSize asize, oarsize,outfsize
COMMON /ERRUNIT/ERRSTOP
INTEGER*8 i,j,k
if(nrhs .gt. 3) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2016:nInput',
& 'Two inputs required.')
elseif(nlhs .gt. 2) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2016:nOutput',
& 'Too many output arguments.')
endif
c get inputs
in_ptr = mxGetPr(prhs(1))
mrows = mxGetM(prhs(1))
ncols = mxGetN(prhs(1))
asize = mrows*ncols
call mxCopyPtrToReal8(in_ptr,ijf,asize)
c copy over jf array
DO 2737 I=1,50
if(ijf(i) .le. 0.0D0) then
jf(i) = .false.
else
jf(i) = .true.
endif
2737 continue
c writes off
jf(38)=.true.
c messages off
jf(34)=.false.
c user inputs
UINPUT = .false.
if( nrhs .gt. 2) then
UINPUT = .true.
uin_ptr = mxGetPr(prhs(3))
umrows = mxGetM(prhs(3))
uncols = mxGetN(prhs(3))
if(umrows.ne.size(uparam)) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2020:numrows',
& 'User parameter array must have size [16,n].')
endif
elseif(.not.all(
c Do not run if requested user inputs are not provided
& JF((/8,9,10,13,14,15,16,17,25,27,32,43,44,45,46/)))) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2020:uimiss',
& 'User Inputs requested but not provided')
endif
nout = 0
if(jf(1)) nout=nout+1
if(jf(2)) nout=nout+3
if(jf(3)) nout=nout+5
if(jf(3).and..not.jf(6)) nout=nout+2
if(.not.jf(24)) nout=nout+1
in_ptr = mxGetPr(prhs(2))
mrows = mxGetM(prhs(2))
ncols = mxGetN(prhs(2))
if(mrows.ne.size(input)) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2020:nmrows',
& 'Date/Location array must have size [9,n].')
endif
if(UINPUT.and.ncols.ne.uncols) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2020:uimismatch',
& 'Date/Location array and user parameter array must match')
endif
asize = 1*mrows
tnmax = 1
nummax=1000
DO 1167 I=1,ncols
call mxCopyPtrToReal8(in_ptr+(8*mrows*(i-1)),input,asize)
heibeg = sngl(input(7))
heiend = sngl(input(8))
heistp = sngl(input(9))
numhei=int(abs(heiend-heibeg)/abs(heistp))+1
if(tnmax.lt.numhei) tnmax=numhei
1167 continue
if(tnmax.gt.nummax) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2020:nummax',
& 'Too many elevations requested.')
endif
nummax = tnmax
c prepare output arrays
outfsize = size(outf);
c plhs(1) = mxCreateDoubleMatrix(size(outf,1)*ncols,nummax,0)
plhs(1) = mxCreateDoubleMatrix(nummax,nout*ncols,0)
outf_ptr = mxGetPr(plhs(1))
if(nlhs .gt. 1) then
plhs(2) = mxCreateDoubleMatrix(size(oar),1*ncols,0)
oar_ptr = mxGetPr(plhs(2))
oarsize=size(oar)
endif
c required by iri_sub
call read_ig_rz
call readapf107
c LOOP OVER ALL INPUTS
DO 1 I=1,ncols
DO 6249 j=1,100
6249 oar(j)=-1.0
if(UINPUT) then
call mxCopyPtrToReal8(uin_ptr+(8*umrows*(i-1)),
& uparam,size(uparam))
C jf(8) =.false. OARR(1)=user input for foF2/MHz or NmF2/m-3
if(.not.jf(8)) then
oar(1)=sngl(uparam(1))
endif
C jf(9) =.false. OARR(2)=user input for hmF2/km or M(3000)F2
if(.not.jf(9)) then
oar(2)=sngl(uparam(2))
endif
C jf(10 )=.false. OARR(15),OARR(16)=user input for Ne(300km),
C Ne(400km)/m-3. Use OARR()=-1 if one of these values is not
C available. If jf(23)=.false. then Ne(300km), Ne(550km)/m-3.
if(.not.jf(10)) then
oar(15)=sngl(uparam(3))
oar(16)=sngl(uparam(4))
endif
C jf(13) =.false. OARR(3)=user input for foF1/MHz or NmF1/m-3
if(.not.jf(13)) then
oar(3)=sngl(uparam(5))
endif
C jf(14) =.false. OARR(4)=user input for hmF1/km
if(.not.jf(14)) then
oar(4)=sngl(uparam(6))
endif
C jf(15) =.false. OARR(5)=user input for foE/MHz or NmE/m-3
if(.not.jf(15)) then
oar(5)=sngl(uparam(7))
endif
C jf(16) =.false. OARR(6)=user input for hmE/km
if(.not.jf(16)) then
oar(6)=sngl(uparam(8))
endif
C jf(17) =.false. OARR(33)=user input for Rz12
if(.not.jf(17)) then
oar(33)=sngl(uparam(9))
endif
C jf(25) =.false. OARR(41)=user input for daily F10.7 index
if(.not.jf(25)) then
oar(41)=sngl(uparam(10))
endif
C jf(27) =.false. OARR(39)=user input for IG12
if(.not.jf(27)) then
oar(39)=sngl(uparam(11))
endif
C jf(32) =.false. OARR(46)=user input for 81-day avg F10.7
if(.not.jf(32)) then
oar(46)=sngl(uparam(12))
endif
C jf(43) =.false. OARR(10)=user input for B0
if(.not.jf(43)) then
oar(10)=sngl(uparam(13))
endif
C jf(44) =.false. OARR(35)=user input for B1
if(.not.jf(44)) then
oar(35)=sngl(uparam(14))
endif
C jf(45) =.false. OARR(89)=user input for HNEA (Ne lower boundary)
if(.not.jf(45)) then
oar(89)=sngl(uparam(15))
endif
C jf(46) =.false. OARR(90)=user input for HNEE (Ne upper boundary)
if(.not.jf(46)) then
oar(90)=sngl(uparam(16))
endif
endif
call mxCopyPtrToReal8(in_ptr+(8*mrows*(i-1)),input,asize)
C Clear error flag
ERRSTOP = -1.0
call IRI_SUB(jf,int(input(1)),sngl(input(2)),sngl(input(3)),
& int(input(4)),int(input(5)),sngl(input(6)),
& sngl(input(7)),
& sngl(input(8)),sngl(input(9)),outf,OAR)
if(ERRSTOP.gt.0) then
call mexErrMsgIdAndTxt ('MATLAB:IRI2020:StopCalled',
& 'IRI_SUB called a stop command. Attempting to exit safely.')
endif
k=0
c output Ne if requested
C OUTF(1,*) ELECTRON DENSITY/M-3
if(jf(1)) then
call mxCopyReal8ToPtr(dble(outf(1,1:nummax)),
& outf_ptr+8*nummax*((i-1)*nout+k),nummax)
k=k+1
endif
c output Ti/Te if requested
C OUTF(2,*) NEUTRAL TEMPERATURE/K
C OUTF(3,*) ION TEMPERATURE/K
C OUTF(4,*) ELECTRON TEMPERATURE/K
if(jf(2)) then
do 43 j=2,4
call mxCopyReal8ToPtr(dble(outf(j,1:nummax)),
& outf_ptr+8*nummax*((i-1)*nout+k),nummax)
k=k+1
43 continue
endif
c output Ni if requested
C OUTF(5,*) O+ ION DENSITY/% or /M-3 if jf(22)=f
C OUTF(6,*) H+ ION DENSITY/% or /M-3 if jf(22)=f
C OUTF(7,*) HE+ ION DENSITY/% or /M-3 if jf(22)=f
C OUTF(8,*) O2+ ION DENSITY/% or /M-3 if jf(22)=f
C OUTF(9,*) NO+ ION DENSITY/% or /M-3 if jf(22)=f
if(jf(3)) then
do 44 j=5,9
call mxCopyReal8ToPtr(dble(outf(j,1:nummax)),
& outf_ptr+8*nummax*((i-1)*nout+k),nummax)
k=k+1
44 continue
endif
c output Ni if requested
C AND, IF JF(6)=.FALSE.:
C OUTF(10,*) CLUSTER IONS DEN/% or /M-3 if jf(22)=f
C OUTF(11,*) N+ ION DENSITY/% or /M-3 if jf(22)=f
if(jf(3).and..not.jf(6)) then
do 45 j=10,11
call mxCopyReal8ToPtr(dble(outf(j,1:nummax)),
& outf_ptr+8*nummax*((i-1)*nout+k),nummax)
k=k+1
45 continue
endif
c output special D-Region values
C if(jf(24) OUTF(14,1:11) standard IRI-Ne for 60,65,..,110km
C =.false.) 12:22) Friedrich (FIRI) model at these heights
C 23:33) standard Danilov (SW=0, WA=0)
C 34:44) for minor Stratospheric Warming (SW=0.5)
C 45:55) for major Stratospheric Warming (SW=1)
C 56:66) weak Winter Anomaly (WA=0.5) conditions
C 67:77) strong Winter Anomaly (WA=1) conditions
if(.not.jf(24)) then
call mxCopyReal8ToPtr(dble(outf(14,1:nummax)),
& outf_ptr+8*nummax*((i-1)*nout+k),nummax)
k=k+1
endif
C OARR(1:100) ADDITIONAL OUTPUT PARAMETERS
if(nlhs .gt. 1) then
call mxCopyReal8ToPtr(dble(oar),
& oar_ptr+(8*(i-1)*oarsize),oarsize)
endif
1 CONTINUE
return
end
三、matlab调用并编译接口函数
IRI2020.m
function Out = IRI2020(lat,lon,time,alt,varargin)
% MATLAB wrapper for the IRI model and MEX function
%
% INPUTS
% latitude(deg), longitude(deg), time(datenum or datetime), altitude(km)
%
% OUTPUTS
% structure with fields
% latitude(deg), longitude(deg), time(datetime), altitude(km),
% electron density (Ne/m^3) and various profile parameters
%
% OPTIONAL INPUTS
%
% 'JF', JF array: use the keyword 'JF' followed by a length 50 numerical
% vector of option switches as defined in irisub.for. The wrapper
% interprets JF<=0 as false, JF>0 as true. By default, the wrapper uses a
% resonable set of options that should work for most purposes. All
% options related to non-Ne outputs or user defined input parameters will
% be ignored.
%
% By default, IRI2016 works in profile mode. Lat and Lon must be the
% same size, and the output will have Ne at every altitude for each
% Lat/Lon pair. Time can either have size 1, and all profiles will be run
% for the same time, or else it can be the same size as lat/lon in which
% case each corresponding [lat,lon,time] will be run
%
% There is no limit on how many altitudes can be given for a profile,
% unlike the base IRI. This function will break calls to IRI_sub up
% if more than 1000 altitudes are requested
%
% 'sat' mode: in this mode, all inputs must be the same size, including
% altitude
%
% 'map' mode: in this mode, all unique values of inputs are expanded in a
% 4d grid. The output Ne is of a size [n time,n lat,n lon,n alt]
ErrorStr = 'MATLAB:IRI2020';
% PATH TO MEX FILE
[RootDir, ~, ~] = fileparts(mfilename('fullpath'));
% PATH TO NRLMSISE2.1 FORTRAN CODE
ModelDir = fullfile(RootDir,'IRI-2020');
ModelDir = fullfile(ModelDir,filesep);
if isempty(ModelDir)
error(ErrorStr,[...
'IRI2016 model directory not configured! Please edit',...
' IRI2016.m to specify the appropriate path'])
elseif ~exist(ModelDir,'dir')
error(ErrorStr,[...
'Invalid model directory. Directory does not exist. Please edit',...
' IRI2016.m to specify the appropriate path'])
end
JF = ones(50,1);
JF([4,5,6,23,30,33,34,35,39,40,47])=-1;
% JF switches to turn off/on (.true./.false.) several options
%
% i .true. .false. standard version
% -----------------------------------------------------------------
% 1 Ne computed Ne not computed t
% 2 Te, Ti computed Te, Ti not computed t
% 3 Ne & Ni computed Ni not computed t
% 4 B0,B1 - Bil-2000 B0,B1 - other models jf(31) false
% 5 foF2 - CCIR foF2 - URSI false
% 6 Ni - DS-1995 & DY-1985 Ni - RBV-2010 & TBT-2015 false
% 7 Ne - Tops: f10.7<188 f10.7 unlimited t
% 8 foF2 from model foF2 or NmF2 - user input t
% 9 hmF2 from model hmF2 or M3000F2 - user input t
% 10 Te - Standard Te - Using Te/Ne correlation t
% 11 Ne - Standard Profile Ne - Lay-function formalism t
% 12 Messages to unit 6 to messages.txt on unit 11 t
% 13 foF1 from model foF1 or NmF1 - user input t
% 14 hmF1 from model hmF1 - user input (only Lay version)t
% 15 foE from model foE or NmE - user input t
% 16 hmE from model hmE - user input t
% 17 Rz12 from file Rz12 - user input t
% 18 IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 t
% 19 F1 probability model only if foF1>0 and not NIGHT t
% 20 standard F1 standard F1 plus L condition t
% (19,20) = (t,t) f1-prob, (t,f) f1-prob-L, (f,t) old F1, (f,f) no F1
% 21 ion drift computed ion drift not computed t
% 22 ion densities in % ion densities in m-3 t
% 23 Te_tops (Bil-1985) Te_topside (TBT-2012) false
% 24 D-region: IRI-1990 FT-2001 and DRS-1995 t
% 25 F107D from APF107.DAT F107D user input (oarr(41)) t
% 26 foF2 storm model no storm updating t
% 27 IG12 from file IG12 - user t
% 28 spread-F probability not computed t
% 29 IRI01-topside new options as def. by JF(30) t
% 30 IRI01-topside corr. NeQuick topside model false
% (29,30) = (t,t) IRIold, (f,t) IRIcor, (f,f) NeQuick, (t,f) IRIcor2
% 31 B0,B1 ABT-2009 B0 Gulyaeva-1987 h0.5 t
% (4,31) = (t,t) Bil-00, (f,t) ABT-09, (f,f) Gul-87, (t,f) not used
% 32 F10.7_81 from file F10.7_81 - user input (oarr(46)) t
% 33 Auroral boundary model on/off true/false false
% 34 Messages on Messages off t
% 35 foE storm model no foE storm updating false
% 36 hmF2 w/out foF2_storm with foF2-storm t
% 37 topside w/out foF2-storm with foF2-storm t
% 38 turn WRITEs off in IRIFLIP turn WRITEs on t
% 39 hmF2 (M3000F2) new models false
% 40 hmF2 AMTB-model Shubin-COSMIC model t
% (39,40) = (t,t) hmF2-old, (f,t) AMTB, (f,f) Shubin, (t,f) not used
% 41 Use COV=F10.7_365 COV=f(IG12) (IRI before Oct 2015) t
% 42 Te with PF10.7 dep. w/o PF10.7 dependance t
% 43 B0 from model B0 user input in OARR(10) t
% 44 B1 from model B1 user input in OARR(35) t
% 45 HNEA=65/80km dya/night HNEA user input in OARR(89) t
% 46 HNEE=2000km HNEE user input in OARR(90) t
% 47 CGM computation on CGM computation off false
% 48 Ti Tru-2021 Bil-1981 t
% ....
% 50
% ------------------------------------------------------------------
mode='default';
Update=false;
Storm=false;
i=1;
while i <= numel(varargin)
if ischar(varargin{i}) || isstring(varargin{i})
if strcmpi(varargin{i},'default')
mode = default;
elseif strcmpi(varargin{i},'sat')
mode = 'sat';
elseif strcmpi(varargin{i},'map')
mode='map';
elseif strcmpi(varargin{i},'update')
Update=true;
elseif strcmpi(varargin{i},'jf')
JF = varargin{i+1};
i=i+1;
else
error(ErrorStr,[...
'Invalid option: ',varargin{i}]);
end
else
error(ErrorStr,[...
'Invalid option.'])
end
i=i+1;
end
if Update % update indices.dat file in iri directory
read_apf(fullfile(ModelDir,'apf107.dat'),true);
read_ig_rz(fullfile(ModelDir,'ig_rz.dat'),true);
end
if any((lat(:)<-90) | (lat(:)>90))
error(ErrorStr,...
'Input latitude out of range [-90,90]')
end
if any((lon(:)<-180) | (lon(:)>360))
error(ErrorStr,...
'Input longitude out of range [-180,360]')
end
try
OldDir = pwd;
cd(ModelDir);
if ~exist('mex_iri2020','file')
warning(ErrorStr,...
'Required MEX file mex_iri2020 does not exist. Attempting to compile.')
mex('-compatibleArrayDims','-DO3','-Dstatic','-DfPIC', ...
'-Dffast-math','-Dmarch=native', ...
'-output',fullfile(RootDir,'mex_iri2020'), ...
fullfile(ModelDir,'irisub.F'), ...
fullfile(ModelDir,'irifun.F'), ...
fullfile(ModelDir,'iritec.F'),...
fullfile(ModelDir,'iridreg.F'), ...
fullfile(ModelDir,'iriflip.F'), ...
fullfile(ModelDir,'cira.F'), ...
fullfile(ModelDir,'igrf.F'),...
fullfile(ModelDir,'rocdrift.F'),...
fullfile(RootDir,'mex_iri.F'))
end
% check to make sure the IRI will have the needed files
if ~all(arrayfun(@(n) exist(sprintf('ccir%02i.asc',n),'file'),11:22))
error(ErrorStr,[...
'CCIR files required by IRI are missing from the model folder',...
' Check to make sure IRI2020.m is configured correctly.'])
end
if ~all(arrayfun(@(n) exist(sprintf('mcsat%02i.dat',n),'file'),11:22))
error(ErrorStr,[...
'mcsat files required by IRI are missing from the model folder',...
' Check to make sure IRI2020.m is configured correctly.'])
end
if ~all(arrayfun(@(n) exist(sprintf('ursi%02i.asc',n),'file'),11:22))
error(ErrorStr,[...
'URSI files required by IRI are missing from the model folder',...
' Check to make sure IRI2020.m is configured correctly.'])
end
if ~exist('ig_rz.dat','file')
error(ErrorStr,[...
'File ig_rz.dat required by IRI is missing from the model folder',...
' Check to make sure IRI2020.m is configured correctly, or', ...
' use the ''update'' keyword to download.'])
end
if ~exist('apf107.dat','file')
error(ErrorStr,[...
'File apf107.dat required by IRI is missing from the model folder',...
' Check to make sure IRI2020.m is configured correctly, or', ...
' use the ''update'' keyword to download.'])
end
time = datetime(datenum(time),'ConvertFrom','datenum');
if strcmpi(mode,'default')
tvec=[numel(lat),numel(lon)];
tvec=tvec==tvec';
if ~all(tvec(:))
error(ErrorStr,[...
'In "default" mode lat and lon must have the same size'])
end
if numel(time) == 1
time = repmat(time, size(lon));
elseif numel(lon) == 1
lon=repmat(lon,size(time));
lat=repmat(lat,size(time));
else
error(ErrorStr,[...
'In "default" mode, the number of input times must be either',...
' equal to one, or the number of lat/lon pairs'])
end
Out.dates = time(:)';
Out.lon = lon(:);
Out.lat = lat(:);
Out.alt = alt(:);
vbeg=[];
vend=[];
vstp=[];
split_alt;
lon = repmat(lon(:)',numel(vstp),1);
lat = repmat(lat(:)',numel(vstp),1);
time = repmat(time(:)',numel(vstp),1);
vbeg = repmat(vbeg(:),1,size(lon,2));
vend = repmat(vend(:),1,size(lon,2));
vstp = repmat(vstp(:),1,size(lon,2));
nout = round((vend-vbeg)./vstp)+1;
nalt = sum(nout(:,1),1);
nummax = max(nout(:));
nper = size(nout,1);
nprof = size(nout,2);
outfIndex = [(1:nalt)']+[(0:(nprof-1)).*nummax.*nper];
OutFFmt = @(X) X(outfIndex);
oarIndex = 1:nper:(nper*nprof);
oarFmt = @(X,io) reshape(X(io,oarIndex),nprof,1);
elseif strcmpi(mode,'sat')
tvec=[numel(lat),numel(lon),numel(alt),numel(time)];
tvec=tvec==tvec';
if ~all(tvec(:))
error(ErrorStr,[...
'In "sat" mode lat, lon, alt, time must all have the same size'])
end
Out.dates = time(:)';
Out.lon = lon(:);
Out.lat = lat(:);
Out.alt = alt(:);
vbeg=alt(:);
vend=alt(:);
vstp=alt(:);
OutFFmt = @(X) X(:);
oarFmt = @(X,io) X(io,:)';
elseif strcmpi(mode,'map')
lat = unique(lat);
lon = unique(lon);
alt = unique(alt);
time = unique(time);
vbeg=[];
vend=[];
vstp=[];
split_alt;
nlon = numel(lon);
nlat=numel(lat);
nalt=numel(alt);
ntime=numel(time);
Out.lat = lat(:);
Out.lon = lon(:);
Out.alt = alt(:);
Out.dates = time(:)';
[vbeg,~,~,~] = ndgrid(vbeg(:),lat(:),lon(:),time(:));
[vend,~,~,~] = ndgrid(vend(:),lat(:),lon(:),time(:));
[vstp,lat,lon,time] = ndgrid(vstp(:),lat(:),lon(:),time(:));
nout = ceil((vend-vbeg)./vstp)+1;
nalt = sum(nout(:,1),1);
nummax = max(nout(:));
nper = size(nout,1);
nprof = nlon*nlat*ntime;
outfIndex = [(1:nalt)']+[(0:(nprof-1)).*nummax.*nper];
OutFFmt = @(X) permute(reshape(X(outfIndex),nalt,nlat,nlon,ntime),[4,2,3,1]);
oarIndex = 1:nper:(nper*nprof);
oarFmt = @(X,io) permute(reshape(X(io,oarIndex),nlat,nlon,ntime),[3,1,2]);
end
noutf = 0;
if JF(1) > 0; noutf = noutf+1; end
if JF(2) > 0; noutf = noutf+3; end
if JF(3) > 0; noutf = noutf+5; end
if JF(3) > 0 && JF(6) < 0; noutf = noutf+2; end
if JF(24) < 0; noutf = noutf+1; end
% no user inputs
JF([8,9,10,13,14,15,16,17,25,27,32,43,44,45,46]) = 1;
% no outputs but our outputs
JF([34,38]) = -1;
UT = 24.*mod(datenum(time(:)),1);
mmdd = 100.*time(:).Month + time(:).Day;
yr = time(:).Year;
ind = [zeros(size(lat(:))),lat(:),lon(:),yr(:),mmdd(:),UT(:)+25,vbeg(:),vend(:),vstp(:)];
[outf,oar] = mex_iri2020(double(JF),double(ind'));
k = 1;
if JF(1) > 0 % output electron density
Out.dens = OutFFmt(outf(:,k:noutf:end));
k = k+1;
end
if JF(2) > 0% output temperatures
fn = {'Tn','Ti','Te'};
for i = 1:numel(fn)
Out.(fn{i}) = OutFFmt(outf(:,k:noutf:end));
k = k+1;
end
end
if JF(3) > 0% output ion densities
fn = {'O_p','H_p','He_p','O2_p','NO_p'};
for i = 1:numel(fn)
Out.(fn{i}) = OutFFmt(outf(:,k:noutf:end));
k = k+1;
end
if JF(6) < 0% output other ions
fn = {'Cluster_p','N_p'};
for i = 1:numel(fn)
Out.(fn{i}) = OutFFmt(outf(:,k:noutf:end));
k = k+1;
end
end
end
iid=[ 1, 2, 3, 4, 5, 6, 10, 35, ...
23, 25, 36, 33, 39, 41, 46, 83];
iin={'NmF2','hmF2','NmF1','hmF1','NmE','hmE','B0','B1', ...
'solzen','dip','M3000F2','Rz12','IG12','F107','F107_81','Kp'};
for ji = 1:numel(iid)
Out.(iin{ji}) = oarFmt(oar,iid(ji));
Out.(iin{ji})(Out.(iin{ji})==-1)=nan;
end
Out.mlat = oarFmt(oar,49);
Out.mlon = oarFmt(oar,50);
Out.mlt = oarFmt(oar,54);
% OARR(1:100) ADDITIONAL OUTPUT PARAMETERS
%
% #OARR(1) = NMF2/M-3 #OARR(2) = HMF2/KM
% #OARR(3) = NMF1/M-3 #OARR(4) = HMF1/KM
% #OARR(5) = NME/M-3 #OARR(6) = HME/KM
% OARR(7) = NMD/M-3 OARR(8) = HMD/KM
% OARR(9) = HHALF/KM #OARR(10) = B0/KM
% OARR(11) =VALLEY-BASE/M-3 OARR(12) = VALLEY-TOP/KM
% OARR(13) = TE-PEAK/K OARR(14) = TE-PEAK HEIGHT/KM
% #OARR(15) = TE-MOD(300KM) #OARR(16) = TE-MOD(400KM)/K
% OARR(17) = TE-MOD(600KM) OARR(18) = TE-MOD(1400KM)/K
% OARR(19) = TE-MOD(3000KM) OARR(20) = TE(120KM)=TN=TI/K
% OARR(21) = TI-MOD(430KM) OARR(22) = X/KM, WHERE TE=TI
% OARR(23) = SOL ZENITH ANG/DEG OARR(24) = SUN DECLINATION/DEG
% OARR(25) = DIP/deg OARR(26) = DIP LATITUDE/deg
% OARR(27) = MODIFIED DIP LAT. OARR(28) = Geographic latitude
% OARR(29) = sunrise/dec. hours OARR(30) = sunset/dec. hours
% OARR(31) = ISEASON (1=spring) OARR(32) = Geographic longitude
% #OARR(33) = Rz12 OARR(34) = Covington Index
% #OARR(35) = B1 OARR(36) = M(3000)F2
% $OARR(37) = TEC/m-2 $OARR(38) = TEC_top/TEC*100.
% #OARR(39) = gind (IG12) OARR(40) = F1 probability
% #OARR(41) = F10.7 daily OARR(42) = c1 (F1 shape)
% OARR(43) = daynr OARR(44) = equatorial vertical
% OARR(45) = foF2_storm/foF2_quiet ion drift in m/s
% #OARR(46) = F10.7_81 OARR(47) = foE_storm/foE_quiet
% OARR(48) = spread-F probability
% OARR(49) = Geomag. latitude OARR(50) = Geomag. longitude
% OARR(51) = ap at current time OARR(52) = daily ap
% OARR(53) = invdip/degree OARR(54) = MLT-Te
% OARR(55) = CGM-latitude OARR(56) = CGM-longitude
% OARR(57) = CGM-MLT OARR(58) = CGM lat eq. aurl bodry
% OARR(59) = CGM-lati(MLT=0) OARR(60) = CGM-lati for MLT=1
% OARR(61) = CGM-lati(MLT=2) OARR(62) = CGM-lati for MLT=3
% OARR(63) = CGM-lati(MLT=4) OARR(64) = CGM-lati for MLT=5
% OARR(65) = CGM-lati(MLT=6) OARR(66) = CGM-lati for MLT=7
% OARR(67) = CGM-lati(MLT=8) OARR(68) = CGM-lati for MLT=9
% OARR(69) = CGM-lati(MLT=10) OARR(70) = CGM-lati for MLT=11
% OARR(71) = CGM-lati(MLT=12) OARR(72) = CGM-lati for MLT=13
% OARR(73) = CGM-lati(MLT=14) OARR(74) = CGM-lati for MLT=15
% OARR(75) = CGM-lati(MLT=16) OARR(76) = CGM-lati for MLT=17
% OARR(77) = CGM-lati(MLT=18) OARR(78) = CGM-lati for MLT=19
% OARR(79) = CGM-lati(MLT=20) OARR(80) = CGM-lati for MLT=21
% OARR(81) = CGM-lati(MLT=22) OARR(82) = CGM-lati for MLT=23
% OARR(83) = Kp at current time OARR(84) = magnetic declination
% OARR(85) = L-value OARR(86) = dipole moment
% OARR(87) = SAX300 OARR(88) = SUX300
% #OARR(89) = HNEA #OARR(90) = HNEE
Out.JF = JF;
cd(OldDir);
catch err
cd(OldDir)
rethrow(err)
end
function split_alt
if numel(alt)>1
dalt = diff(alt(:)');
dalt(abs(dalt)<1e-3) = 0;
hind=find(abs([1,diff(dalt),1])>1e-3);
hind(1)=0;
hind([0,diff(hind)]==1)=[];
while any(diff(hind)>1000)
ib = find(diff(hind)>1000,1);
hind=[hind(1:ib),hind(ib)+999,hind((ib+1):end)];
end
vbeg = alt(hind(1:end-1)+1);
vend = alt(hind(2:end));
vstp = dalt(hind(1:end-1)+1);
else
vbeg=alt;
vend=alt;
vstp=1;
end
end
end两个读取太阳地磁指数的函数：
read_apf.m
function apf = read_apf(file, update, remote_url)
%READ_APF read and update an apf107.dat file as needed for the IRI
%
% Syntax
% ------
% APF = READ_APF(FILE)
% APF = READ_APF(FILE,UPDATE)
% APF = READ_APF(FILE,UPDATE,REMOTE_URL)
%
% Description
% ------------
% APF = READ_APF(FILE) will read the solar and magnetic indices given in
% an apf107.dat file as needed for geophysical models like the
% International Reference Ionosphere (IRI). If no path is given, it will
% look in the current folder by default.
%
% APF = READ_APF(FILE,UPDATE) will attempt to update the existing
% apf107.dat file from a remote source. If the file has already been
% updated today, it will not re-download it. If the download does not
% complete successfully, the existing data will not be overwritten.
%
% APF = READ_APF(FILE,UPDATE,REMOTE_URL) will use a remote source other
% than the default
%
% Input Arguments
% ---------------
%
% Name Description Data Type
% ---- -------------------- ---------
% FILE Path to the apf107.dat file string
%
% UPDATE
% Option to update local file boolean
%
% REMOTE_URL
% URL of remote source data file string
%
% See also http://irimodel.org/indices/IRI-Format-indices-files.pdf
if ~exist('file', 'var') || isempty(file)
file = fullfile(pwd,'apf107.dat');
end
if ~exist('update', 'var') || isempty(update)
update = false;
end
if ~exist('remote_url', 'var') || isempty(remote_url)
remote_url = 'https://chain-new.chain-project.net/echaim_downloads/apf107.dat';
end
if update % don't update if file already updated today
if exist(file,'file')
fdat=dir(file);
if fdat.datenum<(now-1)
update = true;
else
update = false;
end
end
end
apf = [];
while update
update = false;
try
data = webread(remote_url);
catch Err
warning('apf107:DownloadError', ...
'While downloading the indices file the following error occured:\n %s', Err.message)
break
end
try
apf = t_read_apf(data');
catch Err
warning('apf107:CorruptedUpdate',...
'While parsing the updated indices file the following error occured:\n %s', Err.message)
break
end
fid = fopen(file,'w');
fprintf(fid,'%s',data);
fclose(fid);
end
if isempty(apf) % did not update file
apf = t_read_apf();
end
function C = t_read_apf(fp)
formatSpec = '%3d%3d%3d%3d%3d%3d%3d%3d%3d%3d%3d%3d%*3d%5.1f%5.1f%5.1f';
formatSize = cellfun(@str2double,regexp(formatSpec,'(?<=\%\**)[\d]*','match'),'UniformOutput',true);
dw=arrayfun(@(n,m) (1:n) + m,formatSize,[0,cumsum(formatSize(1:end-1))],'UniformOutput',false);
if ~exist('fp','var')
% reads apf107.dat
fp = fopen(file,'r');
D = fread(fp,inf,'uint8=>uint8');
fclose(fp);
else
D = fp;
end
D = char(strread(char(D),'%s','delimiter','\n','whitespace',''));
del = repmat(',',size(D,1),1);
D = cell2mat(cellfun(@(n) [D(:,n),del],dw,'UniformOutput',false));
C = textscan(D',formatSpec,'Delimiter',',');
C = cell2struct(C,{'year','month','day','Ap03','Ap06','Ap09','Ap12','Ap15','Ap18','Ap21','Ap24','Ap','F107','F107_81','F107_365'},2);
C.year=double(C.year)+1900+((double(C.year)<50).*100);
C.date = datetime(C.year,double(C.month),double(C.day));
C = rmfield(C,{'year','month','day'});
end
endread_ig_rz.m
function ig_rz = read_ig_rz(file, update, remote_url)
%READ_IG_RZ read and update an ig_rz.dat file as needed for the IRI
%
% Syntax
% ------
% IG_RZ = READ_IG_RZ(FILE)
% IG_RZ = READ_IG_RZ(FILE,UPDATE)
% IG_RZ = READ_IG_RZ(FILE,UPDATE,REMOTE_URL)
%
% Description
% ------------
% IG_RZ = READ_IG_RZ(FILE) will read the solar and magnetic indices given in
% an ig_rz.dat file as needed for geophysical models like the
% International Reference Ionosphere (IRI). If no path is given, it will
% look in the current folder by default.
%
% IG_RZ = READ_IG_RZ(FILE,UPDATE) will attempt to update the existing
% ig_rz.dat file from a remote source. If the file has already been
% updated today, it will not re-download it. If the download does not
% complete successfully, the existing data will not be overwritten.
%
% IG_RZ = READ_IG_RZ(FILE,UPDATE,REMOTE_URL) will use a remote source other
% than the default
%
% Input Arguments
% ---------------
%
% Name Description Data Type
% ---- -------------------- ---------
% FILE Path to the ig_rz.dat file string
%
% UPDATE
% Option to update local file boolean
%
% REMOTE_URL
% URL of remote source data file string
%
% See also http://irimodel.org/indices/IRI-Format-indices-files.pdf
if ~exist('file', 'var') || isempty(file)
file = fullfile(pwd,'ig_rz.dat');
end
if ~exist('update', 'var') || isempty(update)
update = false;
end
if ~exist('remote_url', 'var') || isempty(remote_url)
remote_url = 'https://chain-new.chain-project.net/echaim_downloads/ig_rz.dat';
end
if update % don't update if file already updated today
if exist(file,'file')
fdat=dir(file);
if fdat.datenum<(now-1)
update = true;
else
update = false;
end
end
end
ig_rz = [];
while update
update = false;
try
data = webread(remote_url);
catch Err
warning('ig_rz:DownloadError', ...
'While downloading the indices file the following error occured:\n %s', Err.message)
break
end
try
ig_rz = t_read_ig_rz(data');
catch Err
warning('ig_rz:CorruptedUpdate',...
'While parsing the updated indices file the following error occured:\n %s', Err.message)
break
end
fid = fopen(file,'w');
fprintf(fid,'%s',data);
fclose(fid);
end
if isempty(ig_rz) % did not update file
ig_rz = t_read_ig_rz();
end
function C = t_read_ig_rz(fp)
if ~exist('fp','var')
% reads ig_rz.dat
fp = fopen(file,'r');
D = fread(fp,inf,'uint8=>uint8');
fclose(fp);
else
D = fp;
end
D = char(strread(char(D),'%s','delimiter','\n','whitespace',''));
Dates = strsplit(D(3,:),',');
DateStart = datetime(str2num(Dates{2}),str2num(Dates{1}),1)-days(14);
DateEnd = datetime(str2num(Dates{4}),str2num(Dates{3}),1) + days(28+14);
Weeks = DateStart:days(7):DateEnd;
Months = unique([Weeks.Year(:),Weeks.Month(:)],'rows');
Q = str2num(reshape(D(5:end,:)',1,[]));
Q = Q(:);
C.IG12 = Q(1:size(Months,1));
C.Rz12 = Q(size(Months,1) + (1:size(Months,1)));
C.date = datetime(Months(:,1),Months(:,2),ones(size(Months(:,1))));
end
end四、主函数
% function main
% IRI2020.m INPUTS
% latitude(deg), longitude(deg), time(datenum or datetime), altitude(km)
%
% OUTPUTS
% structure with fields
% latitude(deg), longitude(deg), time(datetime), altitude(km),
% electron density (Ne/m^3) and various profile parameters
clc;clear;
OutData = IRI2020(36,120,[2019 5 1 4 0 0],100:10:1000);
figure;
plot(OutData.dens,OutData.alt,'k-');
set(gca,'XMinorTick','on','YMinorTick','on');
xlabel('IED/(km)','fontsize',8);ylabel('Height/(km)','fontsize',8);
set(gca,'fontsize',8,'fontname','times')用于国际参考电离层模型的 MATLAB MEX 包装器MEX 文件允许用户从 MATLAB 中执行 FORTRAN 代码。
此存储库添加了三个 MATLAB 函数，一个用于与 MEX 代码交互，另外两个自动下载 apf107.dat 和 ig_rz.dat 文件提供
需要的索引到 IRI MSIS。
IRI2020.m 应在首次运行时自动编译 MEX 文件默认情况下，IRI 将配置为 IRISUB.F 中的默认特定 IRI。有关详细信息，请参阅 IRI2020.m。
fortan源代码和系数：
Indices files at http://irimodel.org/indices/:
ig_rz.dat This file(s) contains the solar and ionospheric indices (IG12, Rz12)
for the time period from Jan 1958 onward. The file is updated
quarterly. It is read by subroutine tcon in irifun.for (ASCII).
[This file will be updated at close to quarterly intervals]
apf107.dat This file provides the 3-hour ap magnetic index and F10.7 daily
81-day and 365-day index from 1960 onward (ASCII).
[This file will be updated at close to quarterly intervals]
Daily updates of these two files are available from the ECHAIM website (David
Themens) as described on irimodel.org.
Coefficients files at http://irimodel.org/COMMON_FILES/:
CCIR%%.dat Monthly coefficient files for the CCIR foF2 and M(3000)F2 models
%%=month+10
URSI%%.dat Monthly coefficient files for the URSI foF2 model
%%=month+10
五、运行测试
完美运行，输出电子密度剖面。
六、代码下载
https://gitcode.com/weixin_44392003/IRI2020_matlab/tree/main参考
- Bilitza D , Altadill D , Zhang Y , et al. The International Reference Ionosphere 2012 – a model of international collaboration[J]. Journal of Space Weather & Space Climate, 2014, 4:689-721.
2 . 国际参考电离层官网International Reference Ionosphere
技术共进，成长同行——讯飞AI开发者社区
更多推荐
34
0- 0
已为社区贡献2条内容
所有评论(0)