- Email: [email protected]
Simulation of ionosphere: a numerical technique
Mathematics and Computers in Simulation XXVI North-Holland
(1984) 1:7-121
SIMCJLATION OF IONOSPHERE: A NUMERICAL TECHNIQUE N.S. CHAUHAN Depurtnratt of Astronony
ctnd Spucu Sumccs. Pun~ah Uttrwrsrty,
Patcda
- 14700,‘. ltldia
A numerical technique has been developed to solve the time-dependent equation of continuity c f plasma for the ionosphere. The solution of the equation gives the behaviour of electron density at any time of the day of the year. A computer program for this has bevn developed which has been used to simulate different physical situations of the ionosphere.
1. introduction
Studies in the ionosphere have an important role in the communication of information through the propagation of radio waves. The propagation of a radio signal after its interaction with the ionosphere depends upon the physical and chemical conditions prevailing in the ionosphere. The reflection of a signal occurs at a height where the frequency of the incident signal (Oi) is of the order of the plasma frequency ( or j of the ionosphere, i.e.
Thus a detailed knowledge of electron density (N) in terms of time of the day and height of the ionized layer becomes important for communication between any two points. There are two ways to know the electron density profile of the ionosphere. One is through a continuous monit0rir.g of the ionosphere using some experimental set up; the other is through a theoretical modelling, which is discussed in this paper. The following time dependent equation of continuity of plasma is solved under Indian conditions (longitude = 75”E)
where Q and L are the production and loss rates of ionization. In the divergence term, V is the transport velocity of plasma. The (r, 6, +) system of coordinates is transformed to a new system (a, x, +) of coordinates such that the new system moves with the motion of magnetic lines of force during the daytime as well as during the nighttime (Baxter and Kendall [l]). The plasma equation of continuity is reduced to a parabolic differential equation which is solved for both Northern and Southern hemispheres. A computer program for the above problem was developed by the author on IBM 360/44. The computer program is general enough to give the electron density of the ionosphere at any time during the day and night, on any day of the year and during any solar activity period. For one complete day calculations, it took about 14-15 minutes on IBM 360/44. The technique has also been us .d to simula,e different phenomena of the ionosphere (Chauhan and Gurm [2,3]).
037%4754/84/$3.00
Q 1984, Elsevier Science Publishers t3.V. (North-Holland)
N.S. Chauhan / Simulation of ionosphere
118
The plasma continuity equation relating the changes in electron density per unit time to production, and divergence of plasma flux, is given by
loss
(3) Following Baxter and Kendall (1968), we transform ( r. 8, 9) coordinates to generalised coordinates system so chosen moves with electromagnetic E x B drift. This transformation is ad%antsgeous particularly at low-latitudes because of the geometry of the magnetic field lines. The equation of a magnetic line of force for dipole field is
(a, s, cp). The new coordinate
r = n sin’8
(4)
where a is radial distance of the point of intersection of the Leld line with the piane of geomagnetic equator. The !ower boundary &) of the problem is taken to be 180 km above the surface of the earth. at the level of maximum production. The height Z is measured upward from this level. s is defined through the relation e - 2,‘2 H,, _ e - Zc/l x2
H,,
(5)
= 1 _ e-%,/2f4,
-
where Z = r - 4) and Ze = a - ro. r. and i-I, are chosen to be 180 km and 85 km respectively above the earth’s surface. Thus, at the lower boundary (Z = 0), x = + 1; at the equator (r = a). s = 0. _Xchanges most rapidly with height in the region where the density of ionization is the greatest. This has an advantage in calculations involving large u. .41so numerical integration can be continued to higher altitudes. if so required. The equilibrium electron density is known to be proportional to exp( - Z/2&) at higher altitudes. We choose 2 more convenient dependent N = 1’.G
where > = exp[ - Z,/2H,,]. If $I is the I = 74, whcye T = 1.37 x 104. is the number of seconds in a radian. After substituting the appropriate values, (3) can be reduced to the following form: i3G
ac
~+TW(+),a=A’
aG
c+B _-
r aG a\-
+
C’G + 3’
(8)
and D’ are constants whose values depend uron the input parameters. Ec,uation (8) can be reduced further to that of the following differential equation whose independent :I i:lWts are a’. x’ and +:
Q hcrc A’. B’, C’
(9)
N.S.
Chauhan
where Q’ is defined as the initial condition
/ Simulafron
of rono.where
(a = a’ at cp= &)) of the ordinary differential
equation
da
-=rW(+).
d+
W(+) is the electromagnetic E x B drift velocity at the geomagnetic equator and for a simple case it can be approximated by a simple harmonic function. The solution of (lo), a = a(+) can be obtained numerically in case IV(+) is not in directly integrable form. Equation (9) is a parabolic differential equation and is solved numerically using the Crank-Nicolson method [S]. Equation (9) is first written in finite differences form, replacing the derivatives by the corresponding approximate difference formula at any point (x,, $Q ). Let two adjacent points on the grid be separated from each other in space direction by an increment Ax and in the time direction by an increment A+ (see Fig. 1). The right-hand side of (9) can be written as
.
,$+Bg+cG] The value P = l/2 scheme:
A$+Bg+cG]+D.
+(1-P)
A+1
is chosen here. Writing
the left-hand
- 2G;+’ + G,$’
I
and right-hand
side in a finite difference
Gh+l _Gktl +B
‘+I
(Ad2 q+
(11)
k
2Ax ‘--I
G,!+l-G;-,
-2G;+G;_, +B (Ax)’
2Ax
+ Cclk+’
I
+ CG:
-I- D,” s
(12)
I
Rearrangrng such that all terms containing G* +’ appear on one side, the above equation becomes AA;
+ I,:++11
+ BBJb + ‘GJk’ ’ + CC;
-+I,:,+,’
= Dn,i;
Fig. 1. The matrix of grid points used in the finite difference scheme.
(13)
N. S. Chauharl / Simulation
120
of ionosphere
where A -----2(dx)z
AA:“=
B
CC;+] = _A,_B 2( Ax)’ D,, =-
’
(Ax)’
’
1
c
1 c DD;+’ =2[~11,, + BD,,] + ‘i: G; + l.y +
4Ax ’
q+ 1- 2G; + G;_,
A “_+--&2’ (Ax)’
W”
4Ax’
D,, =
q+ I - G,t, 2Ax ’
Equation (13) is a system of linear algebraic equations and can be rewritten as AA(J)G(J+l)+BB(J)G(J)+CC(J)G(J-l)=DD(J). Put J= 2, 3, 4 . . . . . 50, then AA(2)G(3)
+ BB(2)G(2) + CC(2)Gtl)
= DD(2)
AA(?)G(3)
+ BB(2)G(2) = DD(1) - CC[2)G( 1).
AA(3)G(4)
+ BB(3)G(3) + CC(3)G(2)
= DD(3).
AA(4)G(5)
+ BB(4)G(4) + CC(4)G(3)
= DD(4),
or
AA(50)G(Sl)
+ BB(50)G(SO) + CC(SO)G(49) = DD(S0).
Since G(51) = G(50), therefore [AA
+ BB(50)]G(50)
+ CC(50)G(39)
= DD(50).
The above equations can be put in the following matrix form: BB(2)
AA(2)
0
CC( 3;
BB(3)
AA(3)
0
.. .
9
G(2)
..-
0
d(3) . .
i
0
0
...
0_
0
.. .
0
CC( 49)
BB(49)
AA(49)
0
CC( 50)
(AA( 50) -t BB( 50))
G(b9) L
r,G ( 50)
DD(2) - CC(2)G(l) =
DD(3) . . DD(49) DD( 50)
The above set of equations can be solved using a factorisation method which utilizes the fact that the matrix elelrlsnt corresponding to this set forms a tridiagonal matrix, i.e. only non-zero elements lie on the main diagonal of the matrix. The solution of the above, matrix gives the value of G’s ard using relation (6). the electron density i A:) can be found out.
N. S. Chuuhan / Simulation of ronosphete
121
We solve the continuity equation only for the northern hemisphere, i.e. for x = 1 to .Y= 0. We make 50 steps in this interval, i.e. 3.x = 0.02. Also A+ = n/180 ( = 4 minutes of time) is taken for speedy convergence of the solution. The boundary condition is G= - D/C
at x = I for all time cp.
The initial condition is assumed to be G = 0 at I = I, (starting time). In (9), B has a singular point at x = 0 and therefore aG/ax must be equal to zero. This condition has beeii incorporated in algorithm. Tht: calculations are started at 0600 LT (I = 0) and end twenty-eight hours later for a(+) greater than r,. The solutio,cl G after a few hours becomes independent of the initial conditions chosen. It took about 14-l 5 minutes on an IBM 360/44 computer for one complete day calculations of the electron density.
Refmmes [ 11 12) [3] (41
R.G. Bdxlcr and P.G. Kendall, Pm. Rqwl SW. Set. , I 304 (1968) 171. N.S. Chouhan and H.S. Gurm, Ind. 1. Radio Space Ph.vs. 8 (1979) 106. N S. C’hauhan nnd H.S. Gurm. J. Amm, Tetr. Ph,w. 42 (1980) 265. J. Crank and P. Nicolson. Pm-. Cunzbtidge Phrlos. Sot. 43 (1947) 50.
(1984) 1:7-121
SIMCJLATION OF IONOSPHERE: A NUMERICAL TECHNIQUE N.S. CHAUHAN Depurtnratt of Astronony
ctnd Spucu Sumccs. Pun~ah Uttrwrsrty,
Patcda
- 14700,‘. ltldia
A numerical technique has been developed to solve the time-dependent equation of continuity c f plasma for the ionosphere. The solution of the equation gives the behaviour of electron density at any time of the day of the year. A computer program for this has bevn developed which has been used to simulate different physical situations of the ionosphere.
1. introduction
Studies in the ionosphere have an important role in the communication of information through the propagation of radio waves. The propagation of a radio signal after its interaction with the ionosphere depends upon the physical and chemical conditions prevailing in the ionosphere. The reflection of a signal occurs at a height where the frequency of the incident signal (Oi) is of the order of the plasma frequency ( or j of the ionosphere, i.e.
Thus a detailed knowledge of electron density (N) in terms of time of the day and height of the ionized layer becomes important for communication between any two points. There are two ways to know the electron density profile of the ionosphere. One is through a continuous monit0rir.g of the ionosphere using some experimental set up; the other is through a theoretical modelling, which is discussed in this paper. The following time dependent equation of continuity of plasma is solved under Indian conditions (longitude = 75”E)
where Q and L are the production and loss rates of ionization. In the divergence term, V is the transport velocity of plasma. The (r, 6, +) system of coordinates is transformed to a new system (a, x, +) of coordinates such that the new system moves with the motion of magnetic lines of force during the daytime as well as during the nighttime (Baxter and Kendall [l]). The plasma equation of continuity is reduced to a parabolic differential equation which is solved for both Northern and Southern hemispheres. A computer program for the above problem was developed by the author on IBM 360/44. The computer program is general enough to give the electron density of the ionosphere at any time during the day and night, on any day of the year and during any solar activity period. For one complete day calculations, it took about 14-15 minutes on IBM 360/44. The technique has also been us .d to simula,e different phenomena of the ionosphere (Chauhan and Gurm [2,3]).
037%4754/84/$3.00
Q 1984, Elsevier Science Publishers t3.V. (North-Holland)
N.S. Chauhan / Simulation of ionosphere
118
The plasma continuity equation relating the changes in electron density per unit time to production, and divergence of plasma flux, is given by
loss
(3) Following Baxter and Kendall (1968), we transform ( r. 8, 9) coordinates to generalised coordinates system so chosen moves with electromagnetic E x B drift. This transformation is ad%antsgeous particularly at low-latitudes because of the geometry of the magnetic field lines. The equation of a magnetic line of force for dipole field is
(a, s, cp). The new coordinate
r = n sin’8
(4)
where a is radial distance of the point of intersection of the Leld line with the piane of geomagnetic equator. The !ower boundary &) of the problem is taken to be 180 km above the surface of the earth. at the level of maximum production. The height Z is measured upward from this level. s is defined through the relation e - 2,‘2 H,, _ e - Zc/l x2
H,,
(5)
= 1 _ e-%,/2f4,
-
where Z = r - 4) and Ze = a - ro. r. and i-I, are chosen to be 180 km and 85 km respectively above the earth’s surface. Thus, at the lower boundary (Z = 0), x = + 1; at the equator (r = a). s = 0. _Xchanges most rapidly with height in the region where the density of ionization is the greatest. This has an advantage in calculations involving large u. .41so numerical integration can be continued to higher altitudes. if so required. The equilibrium electron density is known to be proportional to exp( - Z/2&) at higher altitudes. We choose 2 more convenient dependent N = 1’.G
where > = exp[ - Z,/2H,,]. If $I is the I = 74, whcye T = 1.37 x 104. is the number of seconds in a radian. After substituting the appropriate values, (3) can be reduced to the following form: i3G
ac
~+TW(+),a=A’
aG
c+B _-
r aG a\-
+
C’G + 3’
(8)
and D’ are constants whose values depend uron the input parameters. Ec,uation (8) can be reduced further to that of the following differential equation whose independent :I i:lWts are a’. x’ and +:
Q hcrc A’. B’, C’
(9)
N.S.
Chauhan
where Q’ is defined as the initial condition
/ Simulafron
of rono.where
(a = a’ at cp= &)) of the ordinary differential
equation
da
-=rW(+).
d+
W(+) is the electromagnetic E x B drift velocity at the geomagnetic equator and for a simple case it can be approximated by a simple harmonic function. The solution of (lo), a = a(+) can be obtained numerically in case IV(+) is not in directly integrable form. Equation (9) is a parabolic differential equation and is solved numerically using the Crank-Nicolson method [S]. Equation (9) is first written in finite differences form, replacing the derivatives by the corresponding approximate difference formula at any point (x,, $Q ). Let two adjacent points on the grid be separated from each other in space direction by an increment Ax and in the time direction by an increment A+ (see Fig. 1). The right-hand side of (9) can be written as
.
,$+Bg+cG] The value P = l/2 scheme:
A$+Bg+cG]+D.
+(1-P)
A+1
is chosen here. Writing
the left-hand
- 2G;+’ + G,$’
I
and right-hand
side in a finite difference
Gh+l _Gktl +B
‘+I
(Ad2 q+
(11)
k
2Ax ‘--I
G,!+l-G;-,
-2G;+G;_, +B (Ax)’
2Ax
+ Cclk+’
I
+ CG:
-I- D,” s
(12)
I
Rearrangrng such that all terms containing G* +’ appear on one side, the above equation becomes AA;
+ I,:++11
+ BBJb + ‘GJk’ ’ + CC;
-+I,:,+,’
= Dn,i;
Fig. 1. The matrix of grid points used in the finite difference scheme.
(13)
N. S. Chauharl / Simulation
120
of ionosphere
where A -----2(dx)z
AA:“=
B
CC;+] = _A,_B 2( Ax)’ D,, =-
’
(Ax)’
’
1
c
1 c DD;+’ =2[~11,, + BD,,] + ‘i: G; + l.y +
4Ax ’
q+ 1- 2G; + G;_,
A “_+--&2’ (Ax)’
W”
4Ax’
D,, =
q+ I - G,t, 2Ax ’
Equation (13) is a system of linear algebraic equations and can be rewritten as AA(J)G(J+l)+BB(J)G(J)+CC(J)G(J-l)=DD(J). Put J= 2, 3, 4 . . . . . 50, then AA(2)G(3)
+ BB(2)G(2) + CC(2)Gtl)
= DD(2)
AA(?)G(3)
+ BB(2)G(2) = DD(1) - CC[2)G( 1).
AA(3)G(4)
+ BB(3)G(3) + CC(3)G(2)
= DD(3).
AA(4)G(5)
+ BB(4)G(4) + CC(4)G(3)
= DD(4),
or
AA(50)G(Sl)
+ BB(50)G(SO) + CC(SO)G(49) = DD(S0).
Since G(51) = G(50), therefore [AA
+ BB(50)]G(50)
+ CC(50)G(39)
= DD(50).
The above equations can be put in the following matrix form: BB(2)
AA(2)
0
CC( 3;
BB(3)
AA(3)
0
.. .
9
G(2)
..-
0
d(3) . .
i
0
0
...
0_
0
.. .
0
CC( 49)
BB(49)
AA(49)
0
CC( 50)
(AA( 50) -t BB( 50))
G(b9) L
r,G ( 50)
DD(2) - CC(2)G(l) =
DD(3) . . DD(49) DD( 50)
The above set of equations can be solved using a factorisation method which utilizes the fact that the matrix elelrlsnt corresponding to this set forms a tridiagonal matrix, i.e. only non-zero elements lie on the main diagonal of the matrix. The solution of the above, matrix gives the value of G’s ard using relation (6). the electron density i A:) can be found out.
N. S. Chuuhan / Simulation of ronosphete
121
We solve the continuity equation only for the northern hemisphere, i.e. for x = 1 to .Y= 0. We make 50 steps in this interval, i.e. 3.x = 0.02. Also A+ = n/180 ( = 4 minutes of time) is taken for speedy convergence of the solution. The boundary condition is G= - D/C
at x = I for all time cp.
The initial condition is assumed to be G = 0 at I = I, (starting time). In (9), B has a singular point at x = 0 and therefore aG/ax must be equal to zero. This condition has beeii incorporated in algorithm. Tht: calculations are started at 0600 LT (I = 0) and end twenty-eight hours later for a(+) greater than r,. The solutio,cl G after a few hours becomes independent of the initial conditions chosen. It took about 14-l 5 minutes on an IBM 360/44 computer for one complete day calculations of the electron density.
Refmmes [ 11 12) [3] (41
R.G. Bdxlcr and P.G. Kendall, Pm. Rqwl SW. Set. , I 304 (1968) 171. N.S. Chouhan and H.S. Gurm, Ind. 1. Radio Space Ph.vs. 8 (1979) 106. N S. C’hauhan nnd H.S. Gurm. J. Amm, Tetr. Ph,w. 42 (1980) 265. J. Crank and P. Nicolson. Pm-. Cunzbtidge Phrlos. Sot. 43 (1947) 50.