国际参考电离层2020模型Matlab（IRI-2020）

一、IRI模型

国际参考电离层 （IRI） 是一个国际项目 由空间研究委员会 （COSPAR）和国际无线电科学联盟 （URSI）。该组织组成了一个工作组 （ 60 年代后期至今）产生了电离层的经验标准模型，最新版本是IRI-2020。

对于给定的地点、时间和日期，IRI 模型基本参数如下：

输出参数：
电子密度、电子温度、离子温度、 离子组成（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）

高度范围：
电子密度：白天：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