```python
import numpy as np
from sympy import Rational,symbols,exp,lambdify,sqrt,tanh,log,pi
from RadioArray import RadioArray
from ENUFrame import ENU
import astropy.coordinates as ac
import astropy.units as au
import astropy.time as at
from sympy.vector import CoordSysCartesian


def neQuick(h,No=2.2e11,hmax=368.,Ho=50.):
    '''Ne quick model for one layer'''
    res = np.zeros_like(h)
    g = 1./7.
    rfac = 100.
    dh = h - hmax
    g1 = g * dh
    z = dh / (Ho * (1. + rfac * g1 / (rfac * Ho + g1)))
    ee = np.exp(z)
    dial = ee / 1e7 - 1
    sig = 1./(1. + np.exp(-dial))
    res = No*4.0*ee/(1.0 + ee)**2
    #res[ee>1e7] = No*4.0/ee[ee>1e7]
    #res[ee<=1e7] = No*4.0*ee[ee<=1e7]/(1.0 + ee[ee<=1e7])**2
    #res[z > 40] = 0
    return res

def xef1f2(h):
    f2 = neQuick(h,No=2.2e11,hmax=359.,Ho=46.)
    f1 = neQuick(h,No=3e9,hmax=195.,Ho=20.)
    e = neQuick(h,No=3e9,hmax=90.,Ho=10.)
    return f2 + f1 + e


def symbolicChapman(h, nepeak,hmax,width):
    y = (h-hmax)/width
    return nepeak*exp(Rational(1,2)*(Rational(1) - y - exp(-y)))

def symbolicChapman_def1f2(h = None,cond=0):
    if h is None:
        h = symbols('h')
    if cond == 0:
        p = 1
        w = 1
    if cond == 1:
        p = 5
        w = 1.5
    if cond == 2:
        p = 10
        w = 2.
    #f = symbolicChapman(h,Rational(1e10),Rational(450),Rational(250))
    f2 = symbolicChapman(h,Rational(p*2.2e11),Rational(359),Rational(w*46))
    f1 = symbolicChapman(h,Rational(3e9),Rational(195),Rational(20))
    e = symbolicChapman(h,Rational(3e9),Rational(90),Rational(10))
    #d = symbolicChapman(h,Rational(3e9),Rational(80),Rational(50))
    return f2 + f1 + e  + Rational(4.5e8)
    

def ExampleIRI():
    d = np.genfromtxt('exampleIRI.txt',names=True)
    profile = d['ne']
    return d['height'],d['ne']
    
def plotModels():
    '''Plot various models'''
    import pylab as plt
    h = np.linspace(0,3000,10000)
    ### from IRI
    h1,ne1 = ExampleIRI()
    plt.plot(h1,ne1,c='black',label='Ex. IRI')
    #neQuick layers
    ne3 = xef1f2(h1)
    plt.plot(h1,ne3,c='blue',label='neQuick_ef1f2')
    lam = lambdify(symbols('h'),symbolicChapman_def1f2(),'numpy')
    ne4 = lam(h)
    plt.plot(h,ne4,c='green',label='Chapman_def1f2')
    plt.legend(frameon=False)
    plt.xlabel('Height above surface (km)')
    plt.ylabel(r'Electron density $n_e$ (${\rm m}^{-3}$)')
    #plt.yscale('log')
    plt.grid()
    plt.title('Ionosphere Models')
    plt.show()
    
class Model(object):
    '''Symbolic model object. '''
    def makeOrderList(self,paramDict):
        orderList = []
        for key in paramDict.keys():
            orderList.append(key)
        return orderList
    def makeParamDict(self,paramVec,orderList):
        paramDict = {}
        N = np.size(paramVec)
        i = 0
        while i < N:
            paramDict[orderList[i]] = paramVec[i]
            i += 1
        return paramDict
    def makeParamVec(self,paramDict,orderList):
        N = len(orderList)
        paramVec = np.zeros(N)
        i = 0
        while i < N:
            paramVec[i] = paramDict[orderList[i]]
            i += 1
        return paramVec        

class ElectronContentModel(Model):
    def __init__(self,radioArray = None,**kwargs):
        super(ElectronContentModel,self).__init__(**kwargs)
        if radioArray is None:
            radioArray = RadioArray(arrayFile='arrays/lofar.hba.antenna.cfg')
        self.radioArray = radioArray
        self.enu = ENU(location=radioArray.getCenter().earth_location)
    def saveNeFunc(self,neFunc):
        f = tempfile.SpooledTemporaryFile()
        np.savez(f,neFunc=neFunc)
        f.flush()
        return f       
    
class IriModel(ElectronContentModel):
    def __init__(self,**kwargs):
        '''The ionosphere is modelled as the IRI plus a set of solitons
        with time coherence imposed via linear coherence over short time intervals'''
        super(IriModel,self).__init__(**kwargs)
        self.initIriParams = self.createInitIriParams()
        self.iriParamDict = self.initIriParams.copy()
        self.iriOrder = self.makeOrderList(self.initIriParams)
        self.numIriParams = len(self.iriOrder)
        self.iriFunc = self.makeSymbolicIri()
        print("Generated IRI symbolic function with {0} params".format(self.numIriParams))
        
    def setIriParams(self,paramVec):
        '''Set the paramDict for iri from vector'''
        self.iriParamDict = self.makeParamDict(paramVec,self.iriOrder)
        
    def getIriParams(self):
        return self.makeParamVec(self.iriParamDict,self.iriOrder)
        
    def createInitIriParams(self):
        '''Create an initial random param for a soliton'''
        init = {'fhm' : np.random.uniform(low = 400, high = 600),
                'fw' : np.random.uniform(low = 200,high = 300),
                'nefm' : 10**np.random.uniform(low = 9,high=11),
                'f2hm' : np.random.uniform(low = 300, high = 400),
                'f2w' : np.random.uniform(low = 150,high = 200),
                'nef2m' : 10**np.random.uniform(low = 11,high=12),
                'f1hm' : np.random.uniform(low = 100, high = 140),
                'f1w' : np.random.uniform(low = 60,high = 100),
                'nef1m' : 10**np.random.uniform(low = 9,high=11),
                'ehm' : np.random.uniform(low = 80, high = 120),
                'ew' : np.random.uniform(low = 50,high = 80),
                'neem' : 10**np.random.uniform(low = 9,high=10),
                'dhm' : np.random.uniform(low = 50, high = 80),
                'dw' : np.random.uniform(low = 50,high = 60),
                'nedm' : 10**np.random.uniform(low = 8,high=10)}
        init = {'f2hm' : 360,
                'f2w' : 46,
                'nef2m' : 2.2e11,
                'f1hm' : 195,
                'f1w' : 20,
                'nef1m' : 3e9,
                'ehm' : 90,
                'ew' : 10,
                'neem' : 3e9}
        return init
    
    def makeSymbolicIri(self):
        x,y,z = symbols('x y z')
        #R = Rational(6371)
        R = Rational(int(self.radioArray.getCenter().spherical.distance.to(au.km).value))
        r = sqrt(x**(Rational(2))+y**(Rational(2))+z**(Rational(2)))
        h = r - R
        #f = symbolicChapman(h,*symbols('nefm fhm fw'))
        f2 = symbolicChapman(h,*symbols('nef2m f2hm f2w'))
        f1 = symbolicChapman(h,*symbols('nef1m f1hm f1w'))
        e = symbolicChapman(h,*symbols('neem ehm ew'))
        #d = symbolicChapman(h,*symbols('nedm dhm dw'))
        self.iriFunc = f2 + f1 + e + Rational(4.5e8)
        return self.iriFunc
        
    def generateIri(self,paramVec=None):
        '''Sustitute paramDict into symbolic function'''
        if paramVec is not None:
            self.iriParamDict = self.setIriParams(paramVec)
        self.iriModel = self.iriFunc.subs(self.iriParamDict)
        return self.iriModel
    
    def evaluate(self,X,Y,Z):
        iri = lambdify(symbols('x y z'),self.generateIri(),'numpy')
        ne = iri(X.flatten(),Y.flatten(),Z.flatten()).reshape(X.shape)
        return ne
    
    def plotModel(self):
        import pylab as plt
        iri = lambdify(symbols('x y z'),self.generateIri(),'numpy')
        h = np.linspace(6300,9000,1000)
        plt.plot(h,iri(h,0,0))
        plt.xlabel('geocentric radius (km)')
        plt.ylabel('ne m^-3')
        plt.yscale('log')
        plt.grid()
        plt.show()
        

        
def IntegrateChapman(Npeak,hmax,width):
    y0 = -hmax/width
    ymax = 10#(h+3*w-h)/w
    from scipy.special import erf
    return -Npeak*np.sqrt(2*np.exp(1)*np.pi)*(erf(np.exp(-ymax/2)/np.sqrt(2)) - erf(np.exp(-y0/2)/np.sqrt(2)))
    
def IntegrateChapmanNumerical(Npeak,hmax,width,var=1e10):
    h = np.linspace(0,3000,10000)
    lam = lambdify(symbols('h'),symbolicChapman(symbols('h'), Npeak,hmax,width),'numpy')
    ne = lam(h)
    #import pylab as plt
    #plt.plot(h,ne)
    #plt.show()
    from scipy.integrate import simps
    res = []
    for i in range(100):
        res.append(simps(ne+np.random.normal(loc=0,scale=var,size=h.size),h))
    #import pylab as plt
    #plt.hist(res)
    #plt.show()
    return np.mean(res)

def R2upperbound(frequency=100e6,cond=0):
    n_p = 1.24e-2 * frequency**2
    h = np.linspace(0,3000,10000)
    lam = lambdify(symbols('h'),symbolicChapman_def1f2(h = symbols('h'),cond=cond),'numpy')
    if cond == 0:
        var = 1e10
    if cond == 1:
        var = 1e11
    if cond == 2:
        var = 1e12
    ne = lam(h)
    #import pylab as plt
    #plt.plot(h,ne)
    #plt.show()
    from scipy.integrate import simps
    res = []
    alpha = (ne + var)/n_p
    num = alpha**3
    integrand = num/((1-alpha)**(5./2.))
    print("freq:",frequency,"cond:",cond,"alphahat:{0:1.1e}".format(np.mean(alpha)),"sigma_alpha^2:{0:1.1e}".format(np.mean(alpha**2)-np.mean(alpha)**2))
    return "{0:1.1e}".format(simps(integrand,h)*n_p/8./1e16)




def aPrioriModel(h,zenith):
    '''Return a stratified reference ionosphere based on fitted data depending on zenith angle of the sun in degrees'''
    def peakDensity(n0,dn,tau,b,zenith):
        y = zenith/tau
        return n0 + dn*np.exp(-y**2)/(1. + y**(2*b))
    def peakHeight(z0,dz,rho,chi0,zenith):
        return z0 + dz/(1.+np.exp(-(zenith - chi0)/rho))
    def layerDensity(nm,zm,H,z):
        y = (z - zm)/H
        return nm*np.exp(1./2. * (1. - y - np.exp(-y)))
        
    #D layer
    #nm_d = peakDensity(4e8,5.9e8,58.,1300.,zenith)
    y = zenith/58.
    if y < 1:
        nm_d = 4e8 + 5.9e8*np.exp(-y**2)
    else:
        nm_d = 4e8
    zm_d = peakHeight(81.,7.,7.46,100.,zenith)
    H_d = 8.
    n_d = layerDensity(nm_d,zm_d,H_d,h)
    #E layer
    nm_e = peakDensity(1.6e9,1.6e11,87.,8.7,zenith)
    zm_e = 110.
    H_e = 11.
    n_e = layerDensity(nm_e,zm_e,H_e,h)
    #F1 layer
    nm_f1 = peakDensity(2.0e11,9.1e10,54.,13.6,zenith)
    zm_f1 = 185.
    H_f1 = 40.
    n_f1 = layerDensity(nm_f1,zm_f1,H_f1,h)
    #F2 layer
    nm_f2 = peakDensity(7.7e10,4.4e11,111.,4.8,zenith)
    zm_f2 = peakHeight(242.,75.,7.46,96.,zenith)
    H_f2 = 55.
    n_f2 = layerDensity(nm_f2,zm_f2,H_f2,h)

    n = np.atleast_1d(n_d + n_e + n_f1 + n_f2)
    return n
        

def fitExampleIonospheres(folder):
    import glob
    files = glob.glob("{0}/*".format(folder))
    
    # 0 - DOY
    # 1 - Hour
    # 2 - Solar_zenith_angle, degree
    # 3 - Height, km
    # 4 - Electron_density_Ne, m-3
    # 5 - hmF2, km
    # 6 - hmF1, km
    # 7 - hmE, km
    # 8 - hmD, km
    # 9 - NmF2, m-3
    # 10 - NmF1, m-3
    # 11 - NmE, m-3
    # 12 - NmD, m-3
    # 13 - E_valley_width, km
    d = []
    for f in files:
        d.append(np.genfromtxt(f))
    d = np.vstack(d)
    
    mask = d[:,10] == d[:,10]
    doy = d[mask,0]
    hr = d[mask,1]
    chi = d[mask,2]
    nf2 = d[mask,9]
    nf1 = d[mask,10]
    ne = d[mask,11]
    nd = d[mask,12]
    hf2 = d[mask,5]
    hf1 = d[mask,6]
    he = d[mask,7]
    hd = d[mask,8]
    Z = d[mask,3]
    N = d[mask,4]
    
    
    
    from scipy.optimize import curve_fit
    import pylab as plt
    

    def gauss(x, *p):
        hr = x[0,:]
        doy = x[1,:]
        chi = x[2,:]
        A, mu, sigma, C = p
        return A*np.exp(-(hr-mu)**2/(2.*sigma**2)) + C
    
    def peakDensity(chi,A,C,tau,b):
        return np.abs(C) + (np.abs(A)) *np.exp(-(chi/np.abs(tau))**2)/(1 +  (chi/tau)**(2*np.abs(b)))
    
    def peakDensityFit(chi,*p):
        A,C,tau,b = p
        return peakDensity(chi,A,C,tau,b)
    
    def layerDensity(z,nm,zm,H):
        y = (z-zm)/np.abs(H)
        return np.abs(nm) * np.exp(0.5*(1-y-np.exp(-y)))
    

    def peakHeightFit(pred,*p):
        #doy = pred[0,:]
        #hr = pred[1,:]
        chi = pred[:]
        C,D,ACHI,BCHI = p
        return C + (D)/(1 + np.exp(-(chi-ACHI)/BCHI)) #np.exp((doy - ADOY)**2/(2*BDOY**2)) * np.arctan(BCHI*(chi-ACHI)) * (AHR + BHR*hr)
    
    
    
    p0=(1e8,1e8,45,4)
    coeffnd, var_matrix = curve_fit(peakDensityFit, chi, nd, p0=p0)
    print("nd:",coeffnd)
    plt.scatter(chi,nd,c='blue')
    plt.scatter(chi,peakDensityFit(chi,*coeffnd),c='red')
    plt.show()
    p0=(1e9,1e9,45,4)
    coeffne, var_matrix = curve_fit(peakDensityFit, chi, ne, p0=p0)
    print("ne:",coeffne)
    plt.scatter(chi,ne,c='blue')
    plt.scatter(chi,peakDensityFit(chi,*coeffne),c='red')
    plt.show()
    p0=(1e12,1e11,45,4)
    coeffnf2, var_matrix = curve_fit(peakDensityFit, chi, nf2, p0=p0)
    print("nf2:",coeffnf2)
    plt.scatter(chi,nf2,c='blue')
    plt.scatter(chi,peakDensityFit(chi,*coeffnf2),c='red')
    plt.show()
    p0=(1e12,1e11,45,4)
    mask = hf1 != -1
    coeffnf1, var_matrix = curve_fit(peakDensityFit, chi[mask], nf1[mask], p0=p0)
    print("nf1:",coeffnf1)
    plt.scatter(chi[mask],nf1[mask],c='blue')
    plt.scatter(chi,peakDensityFit(chi,*coeffnf1),c='red')
    plt.show()

    pred = chi
    p0 = (np.min(hd),np.max(hd),100,1)
    plt.scatter(chi,hd)
    plt.scatter(chi,peakHeightFit(pred,*p0),color='red')
    plt.show()
    
    
    coeffhd, var_matrix = curve_fit(peakHeightFit, pred, hd, p0=p0)
    print(coeffhd)
    plt.scatter(doy,hd)
    plt.scatter(doy,peakHeightFit(pred,*coeffhd),color='red')
    plt.show()
    
    p0 = (np.min(he),np.max(he),100,1)
    coeffhe, var_matrix = curve_fit(peakHeightFit, pred, he, p0=p0)
    print(coeffhe)
    plt.scatter(doy,he)
    plt.scatter(doy,peakHeightFit(pred,*coeffhe),color='red')
    plt.show()
    
    p0 = (np.min(hf2),np.max(hf1),100,1)
    coeffhf2, var_matrix = curve_fit(peakHeightFit, pred, hf2, p0=p0)
    print(coeffhf2)
    plt.scatter(doy,hf2)
    plt.scatter(doy,peakHeightFit(pred,*coeffhf2),color='red')
    plt.show()
    
    mask = hf1 != -1
    pred = chi[mask]
    p0 = (np.min(hf1[mask]),np.max(hf1[mask]),100,1)
    coeffhf1, var_matrix = curve_fit(peakHeightFit, pred, hf1[mask], p0=p0)
    print(coeffhf1)
    plt.scatter(doy[mask],hf1[mask])
    plt.scatter(doy[mask],peakHeightFit(pred,*coeffhf1),color='red')
    plt.show()
    

    def DEFDensity(pred,*p):
        chi = pred[0,:]
        z = pred[1,:]
        Hd,He,Hf2,Hf1 = p
        nmd = peakDensityFit(chi,*coeffnd)
        nme = peakDensityFit(chi,*coeffne)
        nmf2 = peakDensityFit(chi,*coeffnf2)
        nmf1 = peakDensityFit(chi,*coeffnf1)
        zmd = peakHeightFit(chi,*coeffhd)
        zme = peakHeightFit(chi,*coeffhe)
        zmf2 = peakHeightFit(chi,*coeffhf2)
        zmf1 = peakHeightFit(chi,*coeffhf1)
        
        return layerDensity(z,nmd,zmd,Hd)+layerDensity(z,nme,zme,He)+layerDensity(z,nmf2,zmf2,Hf2)+layerDensity(z,nmf1,zmf1,Hf1)
        
    pred = np.vstack([chi,Z])
    # p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
    #Hd,He,Hf,Ad,Cd,taud,bd,Ae,Ce,taue,be,Af,Cf,tauf,bf
    p0 = [  8,11,55,40]
    
    coeffH, var_matrix = curve_fit(DEFDensity, pred, N, p0=p0)
    print("H:",coeffH)
    plt.scatter(chi,N)
    plt.scatter(chi,DEFDensity(pred,*coeffH),color='red')
    plt.show()
    plt.scatter(Z,N)
    plt.scatter(Z,DEFDensity(pred,*coeffH),color='red')
    plt.show()
    
if __name__=='__main__':
    h = np.linspace(0,2000,1000)
    import pylab as plt
    plt.plot(h,aPrioriModel(h,0))
    plt.plot(h,aPrioriModel(h,20))
    plt.plot(h,aPrioriModel(h,45))
    plt.plot(h,aPrioriModel(h,90))
    plt.show()
    #fitExampleIonospheres("IriData")
    #print( IntegrateChapman(2.2e12,350,50))
    #print( IntegrateChapmanNumerical(2.2e12,350,50))
    #print(R2upperbound(frequency=25e6,cond=0))
    #print( IntegrateChapman(2.2e12,350,100))
    
    #plotModels()
    #from ENUFrame import *
    #iri = IriModel()
    #print(iri.iriFunc)
    #print (iri.generateIri())
    #xvec = [0]
    #yvec = [0]
    #zvec = np.linspace(0,3000,10000)
    #X,Y,Z = np.meshgrid(xvec,yvec,zvec,indexing='ij')
    #points = ac.SkyCoord(X.flatten()*au.km,Y.flatten()*au.km,Z.flatten()*au.km,frame=iri.enu).transform_to('itrs').cartesian.xyz.to(au.km).value
    
    #for x,y,z in zip(X.flatten(),Y.flatten(),Z.flatten()):
    #    points.append(ac.SkyCoord(x*au.km,y*au.km,y*au.km,frame=iri.enu).transform_to('itrs').cartesian.xyz.to(au.km).value)
    #points = np.array(points)
    #xvec = np.linspace(np.min(points[0,:]),np.max(points[0,:]),len(xvec))
    #yvec = np.linspace(np.min(points[1,:]),np.max(points[1,:]),len(yvec))
    #zvec = np.linspace(np.min(points[2,:]),np.max(points[2,:]),len(zvec))
    #X,Y,Z = np.meshgrid(xvec,yvec,zvec)
    #X = points[0,:].reshape(X.shape)
    #Y = points[1,:].reshape(Y.shape)
    #Z = points[2,:].reshape(Z.shape)
    #X,Y,Z=np.meshgrid(xvec,yvec,zvec,indexing='ij')
    #ne = iri.evaluate(points[0,:],points[1,:],points[2,:])
    #print (ne)
    #import pylab as plt
    #plt.plot(np.linspace(0,3000,10000),ne)
    #plt.show()
    
    
```

a2fcc329ad92f8f0b3bcffa0e47c7c81b20aa114
ipynb
Jupyter Notebook
src/ionotomo/notebooks/IRI.ipynb
Joshuaalbert/IonoTomo
9f50fbac698d43a824dd098d76dce93504c7b879
Apache-2.0
2017-06-22T08:47:07.000Z
2021-07-01T12:33:02.000Z
src/ionotomo/notebooks/IRI.ipynb
Joshuaalbert/IonoTomo
9f50fbac698d43a824dd098d76dce93504c7b879
Apache-2.0
2019-04-03T15:21:19.000Z
2019-04-03T15:48:31.000Z
src/ionotomo/notebooks/IRI.ipynb
Joshuaalbert/IonoTomo
9f50fbac698d43a824dd098d76dce93504c7b879
Apache-2.0
2020-03-01T16:20:00.000Z
2020-07-07T15:09:02.000Z
Qwen/Qwen-72B
1. YES
2. YES