'''Krishan Khurana's magnbetic field model
It is a line-by-line translate of the IDL code in
http://lasp.colorado.edu/mop/resources/code/#khuranamodel
If you publish a paper using the code, it is a good idea
to cite a paper of Khurana & Schwarzl (2005).
Ported to python by Yoshifumi Futaana, 2014.
.. warning::
It is always recommended to use :mod:`irfpy.jupiter.kk2009`.
This :mod:`irfpy.juptier.kk_2009_port` contains a bug,
The keep of this module is kept for the compatibility
with the original IDL code.
**USAGE**:
The :meth:`kk_2009` returns the magnetic field in the system III coordinates.
>>> results = kk_2009(ctimer(2000, 6, 1, 10, 32, 17), 30.00, np.deg2rad(86.88), np.deg2rad(157.94))
>>> print('{r[0]:.4f} {r[1]:.4f} {r[2]:.4f}'.format(r=results[0]))
28.9079 5.0641 -3.5039
>>> print('{:.4f}'.format(results[1]))
-4.1292
>>> print(results[2])
1
.. warning::
Do not use "raw" ctime value.
Call always :func:`ctimer` function where global variables
are set for, likely, leapsecond call...
YF does not like this feature by obvious reasons, but it
is how the original code is implemented.
.. warning::
It looks that the original IDL code has a serious error.
Try the following IDL code in a single session.
::
.compile kk_2009
; Conversion of coordinate is tried
ctime = ctimer(1994, 8, 21, 3, 0, 0)
vecin = [1., 0., 0.]
vecou = [0., 0., 0.]
jrot, 's3c', 'jss', vecin, vecou, ctime
print, 'First answer', vecou
; First answer -0.87623097 -0.48189136 0.0000000
; Here, different time is examined.
ctime = ctimer(2000, 6, 1, 10, 32, 17.0)
kk_2009, ctime, 30.0, 86.88 * 3.1416/ 180., 157.94 * 3.1416/ 180., Brm, Btm, Bpm, ianswer
; Conversion of coordinate exactly same as the first one tried.
ctime = ctimer(1994, 8, 21, 3, 0, 0)
vecin = [1., 0., 0.]
vecou = [0., 0., 0.]
jrot, 's3c', 'jss', vecin, vecou, ctime
print, 'Second answer', vecou
; Second answer -0.92678450 -0.37559352 0.0000000
; The answers differ!!
The module, :mod:`kk_2009_port` is kept the bug as it has been.
This is rather to keep the compatibility (for mainly the purpose of debugging/validating the translation) to IDL code.
On the other hand, YF fixed the bug for the module, :mod:`kk_2009`.
I eliminate using global variables (system variables in IDL original code).
This looks to have been avoided. (Source code lines around 1030).
In practice, using :mod:`kk_2009` is highly recommended.
If you want "bug" version, you simply name kk_2009_port to kk_2009:
.. code-block:: python
> from irfpy.jupiter import kk_2009_port as kk_2009
'''
import logging
logging.basicConfig()
_logger = logging.getLogger('kk2009')
_logger.setLevel(logging.DEBUG)
import datetime
import numpy as np
from scipy import special
### Global variables emulate DEFSYSV in IDL code.
_stheta = None
_sphi = None
_phase = None
_ctime = None
_year = None
_month = None
[docs]def print_global():
print('---')
print('_stheta = {}'.format(_stheta))
print('_sphi = {}'.format(_sphi))
print('_phase = {}'.format(_phase))
print('_ctime = {}'.format(_ctime))
print('_year = {}'.format(_year))
print('_month = {}'.format(_month))
print('---')
[docs]def kk_2009(ctime, R, theta, phi):
''' KK model
:param ctime: J2000 time in seconds (can by calculated by the
ctimer subroutine located in this code)
:param R:
:param theta:
:param phi: Right-handed system III coordinates.
Theta phi in radians.
:returns: A tuple, ((Br, Bt, Bp), ZNS3, mpCheck).
(Br, Bt, Bp) in System III spherical coordinats,
ZNS3 for the distance to the current sheet from the
jovigraphic equator in system III coordinates,
and mpCheck if in magnetopause(=1) or not(=0).
Usage:
First you have to calculate the "ctime", J2000 time in second.
>>> ctime = ctimer(2000, 6, 1, 10, 32, 17.0)
>>> print(ctime)
1086085937.0
Then, call this function
>>> R, theta, phi = 30.00, np.deg2rad(86.88), np.deg2rad(157.94)
>>> Bvec, icheck = kk_2009(ctime, R, theta, phi) # doctest: +SKIP
.. note::
Original comments::
#### ;--------------------------------------------------------------------------------------
#### ;---------------Krishan Khurana's Jupiter magnetic field model for IDL-----------------
#### ;--------------------------------------------------------------------------------------
#### ; NAME:
#### ; kk_2009
#### ;
#### ; PURPOSE:
#### ; Given right-handed system III spherical coordinates, this program
#### ; outputs the strength of jupiter's magnetic field in spherical coordinates
#### ; at the specified coordinates.
#### ;
#### ; DESCRIPTION:
#### ; Calculates the magnetic field of a tilted and warped shielded magnetosphere
#### ; using general deformation technique.
#### ;
#### ; CALLING SEQUENCE:
#### ; kk_2009, ctime, R, theta, phi, Brm, Btm, Bpm, mpCheck
#### ;
#### ; INPUTS:
#### ; ctime - J2000 time in seconds (can by calculated by the
#### ; ctimer subroutine located in this code)
#### ; R, theta, phi - right-handed system III coordinates
#### ;
#### ; OUTPUTS:
#### ; Brm, Btm, Bpm - r, theta, phi components of the magnetic field
#### ; in nanoteslas [nT]
#### ; mpCheck - equals 1 if inside the magnetopause
#### ; - equals 0 if outside the magnetopause
#### ;
#### ; note: as a procedure the output to kk_2009 can easily be changed
#### ; to include more than just the components to the magnetic field
#### ; - for example, the following can also be included:
#### ;
#### ; XJSO, YJSO, ZJSO - x, y, z components of the magnetic field in the
#### ;
#### ; juptier-sun-orbital coordinate system (see JROT subroutine),
#### ; ZNS3 - distance to the current sheet from the
#### ; jovigraphic equator in system III coordinates
#### ;
#### ;
#### ; MODIFICATION HISTORY:
#### ; 2003, written by Krishan Khurana in Fortran 77/90
#### ; 2008, ctimer function added by Mariel Desroche
#### ; (converts year, month, day, hour and second to J2000 seconds)
#### ; 2009, converted to IDL by Adam Shinn
#### ; 2009, leapSecond function added by Adam Shinn and Rob Wilson
#### ; (replaces eetime function to easily append for future leap seconds)
#### ;
#### ; WEBSITE:
#### ; This code is available on the website for the Magnetospheres of
#### ; Outer Planets group at the Univeristy of Colorado at Boulder,
#### ; associated with the Laboratory for Atmospheric and Space Physics:
#### ; http://lasp.colorado.edu/MOP/resources/
#### ;--------------------------------------------------------------------------------------
'''
#### pro kk_2009, ctime, R, theta, phi, Brm, Btm, Bpm, mpCheck; XJSO, YJSO, ZJSO, ZNS3
#### vecin = dblarr(3)
#### vecou = dblarr(3)
vecin = np.zeros([3])
vecou = np.zeros([3])
#### ; constants
#### parmod = dblarr(10)
#### M = 8
#### M1 = 64
#### M2 = 16
#### Mode = 7
#### pi = 3.14159265358979d
#### radian = pi/180d
PARMOD = np.zeros([10])
M = 8
M1 = 64
M2 = 16
Mode = 7
pi = np.pi
radian = pi/180.
#### rho = R*sin(theta)
#### xs3 = rho*cos(phi)
#### ys3 = rho*sin(phi)
#### zs3 = R*cos(theta)
rho = R * np.sin(theta)
xs3 = rho * np.cos(phi)
ys3 = rho * np.sin(phi)
zs3 = R * np.cos(theta)
#### JSun, ctime, stheta, sphi, phase ; given ctime, stheta, sphi, and phase are returned
stheta, sphi, phase = JSun(ctime)
#### defsysv, "!stheta", stheta
#### defsysv, "!sphi", sphi
#### defsysv, "!phase", phase
## I do not understand why DEFSYSV is needed here, but emulated by global variables...
global _stheta
_stheta = stheta
global _sphi
_sphi = sphi
global _phase
_phase = phase
#### ; Now define the nine derivatives
#### ; These are calculated at the original location
#### ; First get the unit normals of the zp axis
#### csheet_N, RNx, RNy, RNz, rho, phi, zs3, stheta, ctime
RNx, RNy, RNz = csheet_N(rho, phi, zs3, stheta, ctime)
####
#### ; The z axis of the mapped location is given by
#### zpx = RNx
#### zpy = RNy
#### zpz = RNz
#### zp = xs3*zpx + ys3*zpy + zs3*zpz
zpx = RNx
zpy = RNy
zpz = RNz
zp = xs3*zpx + ys3*zpy + zs3*zpz
#### ; The y axis of the mapped location is given by
#### ypx = zpy
#### ypy = -zpx
#### ypz = 0d
#### yp1 = sqrt(ypx*ypx + ypy*ypy + ypz*ypz)
#### ypx = ypx/yp1
#### ypy = ypy/yp1
#### ypz = ypz/yp1
#### yp = xs3*ypx + ys3*ypy + zs3*ypz
ypx = zpy
ypy = -zpx
ypz = 0.
yp1 = np.sqrt(ypx*ypx + ypy*ypy + ypz*ypz)
ypx = ypx/yp1
ypy = ypy/yp1
ypz = ypz/yp1
yp = xs3*ypx + ys3*ypy + zs3*ypz
####
#### ; The x axis of the mapped location is given by
#### xpx = ypy*zpz - ypz*zpy
#### xpy = ypz*zpx - ypx*zpz
#### xpz = ypx*zpy - ypy*zpx
#### xp = xs3*xpx + ys3*xpy + zs3*xpz
xpx = ypy*zpz - ypz*zpy
xpy = ypz*zpx - ypx*zpz
xpz = ypx*zpy - ypy*zpx
xp = xs3*xpx + ys3*xpy + zs3*xpz
#### ; Now define the nine derivatives
#### dxpdx = xpx
#### dypdx = ypx
#### dzpdx = zpx
#### dxpdy = xpy
#### dypdy = ypy
#### dzpdy = zpy
#### dxpdz = xpz
#### dypdz = ypz
#### dzpdz = zpz
dxpdx = xpx
dypdx = ypx
dzpdx = zpx
dxpdy = xpy
dypdy = ypy
dzpdy = zpy
dxpdz = xpz
dypdz = ypz
dzpdz = zpz
#### R = sqrt(xs3*xs3 + ys3*ys3 + zs3*zs3)
#### Rp = sqrt(xp*xp + yp*yp + zp*zp)
R = np.sqrt(xs3*xs3 + ys3*ys3 + zs3*zs3)
Rp = np.sqrt(xp*xp + yp*yp + zp*zp)
####
#### ; Now calculate the T matrix
#### Txx = (dypdy*dzpdz - dypdz*dzpdy)
#### Txy = (dxpdz*dzpdy - dxpdy*dzpdz)
#### Txz = (dxpdy*dypdz - dxpdz*dypdy)
#### Tyx = (dypdz*dzpdx - dypdx*dzpdz)
#### Tyy = (dxpdx*dzpdz - dxpdz*dzpdx)
#### Tyz = (dxpdz*dypdx - dxpdx*dypdz)
#### Tzx = (dypdx*dzpdy - dypdy*dzpdx)
#### Tzy = (dxpdy*dzpdx - dxpdx*dzpdy)
#### Tzz = (dxpdx*dypdy - dxpdy*dypdx)
Txx = (dypdy*dzpdz - dypdz*dzpdy)
Txy = (dxpdz*dzpdy - dxpdy*dzpdz)
Txz = (dxpdy*dypdz - dxpdz*dypdy)
Tyx = (dypdz*dzpdx - dypdx*dzpdz)
Tyy = (dxpdx*dzpdz - dxpdz*dzpdx)
Tyz = (dxpdz*dypdx - dxpdx*dypdz)
Tzx = (dypdx*dzpdy - dypdy*dzpdx)
Tzy = (dxpdy*dzpdx - dxpdx*dzpdy)
Tzz = (dxpdx*dypdy - dxpdy*dypdx)
#### ; The mapped location is:
#### vecin = [xs3, ys3, zs3]
#### JROT, 's3c', 'jso', vecin, vecou, ctime
#### XJSO = vecou[0]
#### YJSO = vecou[1]
#### ZJSO = vecou[2]
vecin = [xs3, ys3, zs3]
vecou = jrot('s3c', 'jso', vecin, ctime)
XJSO = vecou[0]
YJSO = vecou[1]
ZJSO = vecou[2]
ZNS3 = csheet_struc(rho, phi, XJSO, YJSO, stheta, ctime)
####
#### RLT = atan(YJSO, XJSO)
#### rmap = sqrt(xp*xp + yp*yp)
#### pmap = atan(yp, xp)
#### zmap = zs3 - ZNS3
RLT = np.arctan2(YJSO, XJSO)
rmap = np.sqrt(xp*xp + yp*yp)
pmap = np.arctan2(yp, xp)
zmap = zs3 - ZNS3
####
#### ; Now calculate the field at the mapped location in cylindrical coordinates
#### scol = pi/2d - stheta
scol = pi/2. - stheta
####
#### Using sample.pro, the current values are:
# scol = 1.5163628
# sphi = 2.7565556
# scolOut = <UNDEF>
# sphiOut = <UNDEF>
# xpx, xpy, xpz, ypx, ypy, ypz, zpx, zpy, zpz =
# 0.67232946 -0.73637952 0.075619499 0.73849401 0.67426003
# 0.0000000 -0.050987206 0.055844548 0.99713675
# So far so good. 2014-05-04T11:50:09 (3-significant digits)
#### get_mapped_sunangle, scol, sphi, scolOut, sphiOut, xpx, xpy, xpz, ypx, ypy, ypz, zpx, zpy, zpz
####
scolOut, sphiOut = get_mapped_sunangle(scol, sphi, xpx, xpy, xpz, ypx, ypy, ypz, zpx, zpy, zpz)
#### IDL returns
# scolOut, sphiOut = 0.12268640 -2.6928876
#### Python returns
# So far so good. 2014-05-04T11:58:45
#################
#### Using sample.pro, the current values are:
# PARMOD = 0., 0. ....
# ctime = 1.0860859e+09
# rmap, pmap, zmap =
# 29.774504 -2.6928876 5.7614139
# Brds, Bpds, Bzds = UNDEF
# sphiOut = -2.6928876
#### Using python,
# ctime = 1086085937.0
# rmap, pmap, zmap = (29.774497437565078, -2.6928974483707622, 5.7614152250836614)
# sphiOut = -2.6928974483707622
#### So far so good. 2014-05-04T12:03:25
#### dipole_shield_cyl_S3, PARMOD, ctime, rmap, pmap, zmap, Brds, Bpds, Bzds, sphiOut
(Brds, Bpds, Bzds) = dipole_shield_cyl_S3(PARMOD, ctime, rmap, pmap, zmap, sphiOut)
#### Using sample.pro, the current values are:
# Brds, Bpds, Bzds = 8.5278990 -1.2858312e-06 -14.171122
#### Using python,
# DEBUG:kk2009:Brds=8.527907547774404
# DEBUG:kk2009:Bpds=-1.2850552968401985e-06
# DEBUG:kk2009:Bzds=-14.171125158004727
#### So far so good. 2014-05-04T14:09:02
_logger.debug('Brds={}'.format(Brds))
_logger.debug('Bpds={}'.format(Bpds))
_logger.debug('Bzds={}'.format(Bzds))
#### tail_mag_notilt_all_modes, rmap, pmap, zmap, RLT, Brcs, Bpcs, Bzcs, Mode
Brcs, Bpcs, Bzcs = tail_mag_notilt_all_modes(rmap, pmap, zmap, RLT, Mode)
#### tail_mag_shield_cyl_S3, ctime, rmap, pmap, zmap, Brcss, Bpcss, Bzcss, M, sphiOut, Mode
Brcss, Bpcss, Bzcss = tail_mag_shield_cyl_S3(ctime, rmap, pmap, zmap, M, sphiOut, Mode)
####
### sample.pro
# Brcss, Bpcss, Bzcss = -0.17749937 -0.023122741 -2.8439913
### python
# >>> Brcss, Bpcss, Bzcss
# (-0.17749938914011432, -0.023122815271786182, -2.8439910833035871
# Sofarsogood. 2014-05-04T19:17:32
#### Brmap = Brds + Brcs + Brcss
#### Bpmap = Bpds + Bpcs + Bpcss
#### Bzmap = Bzds + Bzcs + Bzcss
Brmap = Brds + Brcs + Brcss
Bpmap = Bpds + Bpcs + Bpcss
Bzmap = Bzds + Bzcs + Bzcss
####
#### cyl2car_vect, Brmap, Bpmap, Bzmap, BXcarMap, BYcarMap, BZcarMap, pmap
BXcarMap, BYcarMap, BZcarMap = cyl2car_vect(Brmap, Bpmap, Bzmap, pmap)
####
#### cyl2car_vect, Brds, Bpds, Bzds, BXdsMap, BYdsMap, BZdsMap, pmap
BXdsMap, BYdsMap, BZdsMap = cyl2car_vect(Brds, Bpds, Bzds, pmap)
#### cyl2car_vect, Brcs, Bpcs, Bzcs, BXcsMap, BYcsMap, BZcsMap, pmap
BXcsMap, BYcsMap, BZcsMap = cyl2car_vect(Brcs, Bpcs, Bzcs, pmap)
#### cyl2car_vect, Brcss, Bpcss, Bzcss, BXcssMap, BYcssMap, BZcssMap, pmap
BXcssMap, BYcssMap, BZcssMap = cyl2car_vect(Brcss, Bpcss, Bzcss, pmap)
####
#### Bxfinal = Txx*BXcarMap + Txy*BYcarMap + Txz*BZcarMap
#### Byfinal = Tyx*BXcarMap + Tyy*BYcarMap + Tyz*BZcarMap
#### Bzfinal = Tzx*BXcarMap + Tzy*BYcarMap + Tzz*BZcarMap
Bxfinal = Txx*BXcarMap + Txy*BYcarMap + Txz*BZcarMap
Byfinal = Tyx*BXcarMap + Tyy*BYcarMap + Tyz*BZcarMap
Bzfinal = Tzx*BXcarMap + Tzy*BYcarMap + Tzz*BZcarMap
####
#### Bxdsfinal = Txx*BXdsMap + Txy*BYdsMap + Txz*BZdsMap
#### Bydsfinal = Tyx*BXdsMap + Tyy*BYdsMap + Tyz*BZdsMap
#### Bzdsfinal = Tzx*BXdsMap + Tzy*BYdsMap + Tzz*BZdsMap
Bxdsfinal = Txx*BXdsMap + Txy*BYdsMap + Txz*BZdsMap
Bydsfinal = Tyx*BXdsMap + Tyy*BYdsMap + Tyz*BZdsMap
Bzdsfinal = Tzx*BXdsMap + Tzy*BYdsMap + Tzz*BZdsMap
####
#### Bxcsfinal = Txx*BXcsMap + Txy*BYcsMap + Txz*BZcsMap
#### Bycsfinal = Tyx*BXcsMap + Tyy*BYcsMap + Tyz*BZcsMap
#### Bzcsfinal = Tzx*BXcsMap + Tzy*BYcsMap + Tzz*BZcsMap
Bxcsfinal = Txx*BXcsMap + Txy*BYcsMap + Txz*BZcsMap
Bycsfinal = Tyx*BXcsMap + Tyy*BYcsMap + Tyz*BZcsMap
Bzcsfinal = Tzx*BXcsMap + Tzy*BYcsMap + Tzz*BZcsMap
####
#### Bxcssfinal = Txx*BXcssMap + Txy*BYcssMap + Txz*BZcssMap
#### Bycssfinal = Tyx*BXcssMap + Tyy*BYcssMap + Tyz*BZcssMap
#### Bzcssfinal = Tzx*BXcssMap + Tzy*BYcssMap + Tzz*BZcssMap
Bxcssfinal = Txx*BXcssMap + Txy*BYcssMap + Txz*BZcssMap
Bycssfinal = Tyx*BXcssMap + Tyy*BYcssMap + Tyz*BZcssMap
Bzcssfinal = Tzx*BXcssMap + Tzy*BYcssMap + Tzz*BZcssMap
####
#### CAR2SPH_MAG, Bxfinal, Byfinal, Bzfinal, Brm, Btm, Bpm, theta, phi
Brm, Btm, Bpm = CAR2SPH_MAG(Bxfinal, Byfinal, Bzfinal, theta, phi)
####
#### CAR2SPH_MAG, Bxdsfinal, Bydsfinal, Bzdsfinal, Brdsm, Btdsm, Bpdsm, theta, phi
Brdsm, Btdsm, Bpdsm = CAR2SPH_MAG(Bxdsfinal, Bydsfinal, Bzdsfinal, theta, phi)
####
#### CAR2SPH_MAG, Bxcsfinal, Bycsfinal, Bzcsfinal, Brcsm, Btcsm, Bpcsm, theta, phi
Brcsm, Btcsm, Bpcsm = CAR2SPH_MAG(Bxcsfinal, Bycsfinal, Bzcsfinal, theta, phi)
####
#### CAR2SPH_MAG, Bxcssfinal, Bycssfinal, Bzcssfinal, Brcssm, Btcssm, Bpcssm, theta, phi
Brcssm, Btcssm, Bpcssm = CAR2SPH_MAG(Bxcssfinal, Bycssfinal, Bzcssfinal, theta, phi)
####
#### JOVIAN_VIP4_no_dipole, 4, R, theta, phi, BR, BT, BF
BR, BT, BF = JOVIAN_VIP4_no_dipole(4, R, theta, phi)
# IDL: 0.11955383 0.088434597 -0.091125108
# BR, BT, BF=0.11955384397984081 0.08843459277927272 -0.09112511354623398
_logger.debug('BR, BT, BF={} {} {}'.format(BR, BT, BF))
#### Brm = Brm + BR
#### Btm = Btm + BT
#### Bpm = Bpm + BF
Brm = Brm + BR
Btm = Btm + BT
Bpm = Bpm + BF
####
#### Brdsm = Brdsm + BR
#### Btdsm = Btdsm + BT
#### Bpdsm = Bpdsm + BF
Brdsm = Brdsm + BR
Btdsm = Btdsm + BT
Bpdsm = Bpdsm + BF
####
#### ; now calculate the solar wind magnetic field
#### getBimf, ctime, theta, phi, BrS3IMF, BtS3IMF, BpS3IMF
BrS3IMF, BtS3IMF, BpS3IMF = getBimf(ctime, theta, phi)
_logger.debug('BiS3IMF = {} {} {}'.format(BrS3IMF, BtS3IMF, BpS3IMF))
####
#### CheckIfInsideMappedMp, ctime, XS3, YS3, ZS3, ZNS3, mpCheck
mpCheck = CheckIfInsideMappedMp(ctime, xs3, ys3, zs3, ZNS3)
#### ;checkIfInsideMagnetop, ctime, XS3, YS3, ZS3, mpCheck
####
#### ;print, BrS3IMF, BtS3IMF, BpS3IMF
#### if (mpCheck eq 1) then begin
if mpCheck == 1:
#### Brm = Brm + BrS3IMF
#### Btm = Btm + BtS3IMF
#### Bpm = Bpm + BpS3IMF
Brm = Brm + BrS3IMF
Btm = Btm + BtS3IMF
Bpm = Bpm + BpS3IMF
#### endif else begin
else:
#### Brm = BrS3IMF
#### Btm = BtS3IMF
#### Bpm = BpS3IMF
Brm = BrS3IMF
Btm = BtS3IMF
Bpm = BpS3IMF
#### endelse
return (Brm, Btm, Bpm), ZNS3, mpCheck
#### end ; kkk_2008 main program
####
################# Ported.
#### ;--------------------------------------------------------------------------------------
#### ; function doy
#### ;
#### ; day of year calculator
#### ;--------------------------------------------------------------------------------------
#### function doy, iyr, imon, iday
[docs]def doy(iyr, imon, iday):
''' DOY calculator
I like simpler calculation using ``datetime``, but keep it as it has been...
>>> print(doy(1996, 3, 1))
61
>>> print(doy(1999, 3, 1))
60
>>> print(doy(2000, 3, 1))
61
>>> print(doy(2100, 3, 1))
60
'''
#### mon = [31,28,31,30,31,30,31,31,30,31,30,31]
mon = [31,28,31,30,31,30,31,31,30,31,30,31]
#### dy = 0L
dy = 0
####
#### for i = 2, imon do begin
#### dy = dy + mon[i - 2]
for i in range(2, imon + 1):
dy = dy + mon[i - 2]
#### ; Add an extra day for February
#### if (i eq 3) then begin
## !!!!! YF does not remember if || and && in IDL is priority-order (&& first than ||)
## !!!!! or left-to-right order, but if the first case (as of python and C), this logic
## !!!!! could be wrong? **(FIRST || SECOND) && THRID** is correct, **FIRST || SECOND && THRID** is wrong.
#### if ((iyr mod 100) ne 0) || ((iyr mod 400) eq 0) && ((iyr mod 4) eq 0) then dy = dy + 1
#### endif
if i == 3:
if (((iyr % 100) != 0) or ((iyr % 400) == 0)) and ((iyr % 4) == 0):
dy = dy + 1
#### endfor
####
#### dy = dy + iday
dy = dy + iday
#### return, dy
return dy
#### end
####
################# Ported.
#### ;--------------------------------------------------------------------------------------
#### ; function ctimer
#### ;
#### ; The following lines (and others between comment bars below) were added
#### ; by M. Desroche for use with data sets that use DOY instead of mon/day format
#### ; The similiar lines above need to be commented out and these lines commented in
#### ; for use with this data format
#### ;**************************************************************
#### ; function ctimer, iyr, iday, ihr, imin, sec
#### ;**************************************************************
#### ;--------------------------------------------------------------------------------------
#### function ctimer, iyr, imon, iday, ihr, imin, sec
[docs]def ctimer(iyr, imon, iday, ihr, imin, sec):
''' Calculate J2000 time.
If YF understand correctly, this is the epoch from 1966-01-01T00:00:00.
In this case, we can make simpler function using dateutil or time modules.
>>> print(ctimer(1966, 1, 1, 0, 0, 0))
0.0
>>> print(ctimer(1966, 1, 1, 0, 0, 1))
1.0
>>> print(ctimer(1967, 1, 1, 0, 0, 0))
31536000.0
The following test is in sample.pro.
>>> print(ctimer(2000, 6, 1, 10, 32, 17))
1086085937.0
'''
#### ;dayofyear = julday(imon, iday, iyr, 0, 0, 0) - julday(12, 31, iyr-1, 0, 0, 0)
# dayofyear = datetime.datetime(iyr, imon, iday) - datetime.datetime(iyr - 1, 12, 31)
# dayofyear = dayofyear.days
#### ctime = 0d
#### ndays = 0L
ctime = 0.
ndays = 0
#### ; First calculate the number of days from Jan 1, 1966
#### if (iyr gt 1966) then begin
if iyr > 1966:
#### for i = 1966, iyr - 1 do begin
#### ndays = ndays + 365
#### if ((i mod 100) ne 0) || ((i mod 400) eq 0) && ((i mod 4) eq 0) then ndays = ndays + 1
#### endfor
for i in range(1966, iyr):
ndays += 365
if (((i % 100) != 0) or ((i % 400) == 0)) and ((i % 4) == 0):
ndays += 1
#### ; Now add the number of days of the current year
#### ndays = ndays + doy(iyr, imon, iday) - 1
ndays += (doy(iyr, imon, iday) - 1)
#### ;**************************************************************
#### ; ndays = ndays + iday - 1
#### ;**************************************************************
####
#### ctime = ndays*86400d + ihr*3600d + imin*60d + sec
ctime = ndays * 86400. + ihr * 3600. + imin * 60. + sec
#### endif else begin
else:
#### ; Calculate the seconds for the negative years
#### for i = iyr, 1965 do begin
#### ndays = ndays - 365
#### if ((i mod 100) ne 0) || ((i mod 400) eq 0) && ((i mod 4) eq 0) then ndays = ndays - 1
#### endfor
for i in range(iyr, 1966):
ndays = ndays - 365
if (((i % 100) != 0) or ((i % 400) == 0)) and ((i % 4) == 0):
ndays = ndays - 1
#### ; Now subtract the number of days of the current year
#### ndays = ndays + doy(iyr, imon, iday)
ndays = ndays + doy(iyr, imon, iday)
#### ;***************************************************************
#### ; ndays = ndays + iday
#### ;***************************************************************
####
#### ctime = (ndays - 1)*86400d + ihr*3600d + imin*60d + sec
ctime = (ndays - 1) * 86400. + ihr * 3600. + imin * 60. + sec
#### endelse
####
#### ;jultime = julday(imon, iday, iyr, ihr, imin, sec)*24d*3600d
#### ;print, "ctime", ctime
#### ;print, "jultime", jultime
####
#### defsysv, "!ctime", ctime
#### defsysv, "!year", iyr
#### defsysv, "!month", imon
global _ctime
_ctime = ctime
global _year
_year = iyr
global _month
_month = imon
#### return, ctime
return ctime
#### end
####
################# Ported.
#### ;--------------------------------------------------------------------------------------
#### ; function leapSecond(ctime)
#### ; allows ctime to account for leap seconds
#### ;--------------------------------------------------------------------------------------
#### function leapSecond, ctime
[docs]def leapSecond(ctime):
''' If leap second identified, consider it.
Indeed, the implementation uses global variables. Very quick way.
What is ctime here? In the end, the returned is ctime - 34year.
ctime is 2000-01-01?
.. todo::
Look for definition of J2000 epoch!
::
> kk.leapSecond(kk.ctimer(2000, 1, 1, 0, 0, 0))
-43135.81599998474
'''
#### time = !year*100L + !month ; formats year and month as YYYYMM
time = _year * 100 + _month # Year month in YYYYMM format
#### if (time lt 197207) then Tcor = 10d else $ ; leap second jun 30 1972
if time < 197207:
Tcor = 10.
#### if (time lt 197301) then Tcor = 11d else $ ; leap second dec 31 1972
elif time < 197301:
Tcor = 11.
#### if (time lt 197401) then Tcor = 12d else $ ; leap second dec 31 1973
elif time < 197401:
Tcor = 12.
#### if (time lt 197501) then Tcor = 13d else $ ; leap second dec 31 1974
elif time < 197501:
Tcor = 13.
#### if (time lt 197601) then Tcor = 14d else $ ; leap second dec 31 1975
elif time < 197601:
Tcor = 14.
#### if (time lt 197701) then Tcor = 15d else $ ; leap second dec 31 1976
elif time < 197701:
Tcor = 15.
#### if (time lt 197801) then Tcor = 16d else $ ; leap second dec 31 1977
elif time < 197801:
Tcor = 16.
#### if (time lt 197901) then Tcor = 17d else $ ; leap second dec 31 1978
elif time < 197901:
Tcor = 17.
#### if (time lt 198001) then Tcor = 18d else $ ; leap second dec 31 1979
elif time < 198001:
Tcor = 18.
#### if (time lt 198107) then Tcor = 19d else $ ; leap second jun 30 1981
elif time < 198107:
Tcor = 19.
#### if (time lt 198207) then Tcor = 20d else $ ; leap second jun 30 1982
elif time < 198207:
Tcor = 20.
#### if (time lt 198307) then Tcor = 21d else $ ; leap second jun 30 1983
elif time < 198307:
Tcor = 21.
#### if (time lt 198507) then Tcor = 22d else $ ; leap second jun 30 1985
elif time < 198507:
Tcor = 22.
#### if (time lt 198801) then Tcor = 23d else $ ; leap second dec 31 1987
elif time < 198801:
Tcor = 23.
#### if (time lt 199001) then Tcor = 24d else $ ; leap second dec 31 1989
elif time < 199001:
Tcor = 24.
#### if (time lt 199101) then Tcor = 25d else $ ; leap second dec 31 1990
elif time < 199101:
Tcor = 25.
#### if (time lt 199207) then Tcor = 26d else $ ; leap second jun 30 1992
elif time < 199207:
Tcor = 26.
#### if (time lt 199307) then Tcor = 27d else $ ; leap second jun 30 1993
elif time < 199307:
Tcor = 27.
#### if (time lt 199407) then Tcor = 28d else $ ; leap second jun 30 1994
elif time < 199407:
Tcor = 28.
#### if (time lt 199601) then Tcor = 29d else $ ; leap second dec 31 1995
elif time < 199601:
Tcor = 29.
#### if (time lt 199707) then Tcor = 30d else $ ; leap second jun 30 1997
elif time < 199707:
Tcor = 30.
#### if (time lt 199901) then Tcor = 31d else $ ; leap second dec 31 1998
elif time < 199901:
Tcor = 31.
#### if (time lt 200601) then Tcor = 32d else $ ; leap second dec 31 2005
elif time < 200601:
Tcor = 32.
#### if (time lt 201007) then Tcor = 33d else begin ; next potential leap second
elif time < 201007:
Tcor = 33.
else:
logger.warn('Time is out of range for leapSecond function')
raise NotImplementedError('As in the case for IDL, the execution stops here...')
Tcor = 33. # So far, use 33....
#### print, "Time is out of range for leapSecond function"
#### STOP
#### endelse
#### thirty_four_years = 0.1072958367816d10
thirty_four_years = 0.1072958367816e10
#### leaptime = ctime + Tcor - thirty_four_years ; sets time to J2000 epoch
leaptime = ctime + Tcor - thirty_four_years # Set time to J2000 epoch
return leaptime
#### return, leaptime
#### end
####
[docs]def JSun(ctime):
''' (To be written)
:param ctime: J2000 time in seconds (can by calculated by the
ctimer subroutine located in this code)
:returns (stheta, sphi, phase):
stheta, sphi is latitude and longitude (in radians) of the Sun in system III (RH).
phase is Orbital phase angle of Jupiter (in radians)
.. todo::
Carefully address the difference of IDL's mod and python's %.
For IDL, neg mod pos can return neg, but not for python. (Really?)
'''
#### ;--------------------------------------------------------------------------------------
#### ; procedure JSun
#### ;
#### ; INPUT: ctime of the data point
#### ; OUTPUTS: stheta, sphi, latitude and longitude (in radians) of the Sun in system III (RH).
#### ; OUTPUTS: phase, Orbital phase angle of Jupiter (in radians)
#### ; The equations are written in eetime the J200 time convention followed by PDS.
#### ; We first convert ctime to eetime
#### ;--------------------------------------------------------------------------------------
#### pro JSun, ctime, stheta, sphi, phase
####
#### ; Last updated August 12, 2002
#### ; The program first calculates the direction of the sun in non-rotating
#### ; coordinates from equations of the type:
#### ; theta=a1*cos(omegay*t)+a2*sin(omegay*t)+a3*cos(2.*omegay*t)+a4*sin(2.*omegay*t)+a5
#### ; fphi=b1*cos(omegay*t)+b2*sin(omegay*t)+b3*cos(2.*omegay*t)+b4*sin(2.*omegay*t)+b5
#### ; Then we rotate into the System III coordinates
####
#### ; ctime is in double precision. theta and phi are single precision variables
#### ; omega is jupiter's rotation rate
#### ; omegay is jupiter's yearly orbital rate
####
#### ; Initialize variables
#### pi = 3.14159265358979d
#### twopi = 2d*PI
#### radian = PI/180d
#### degree = 180d/PI
#### yrjup = 11.85652502d*86400d*365.25d
#### omega = 870.536d/86400d
#### omegay = 2*PI/yrjup
#### year = 86400d*0.36525d3
#### D360 = 360d
#### etime1 = -8.25767955817479d8
#### three = 3.123d*radian
#### tan3 = 0.054560676d
#### sin3 = 0.054479647d
#### cos3 = 0.99851488d
pi = np.pi
twopi = pi * 2.
radian = pi / 180.
degree = 180. / pi
yrjup = 11.85652502 * 86400 * 365.25
omega = 870.536 / 86400
omegay = 2* pi / yrjup
year = 86400 * 365.25
D360 = 360.
etime1 = -8.25767955817479e8
three = 3.123 * radian
tan3 = 0.054560676
sin3 = 0.054479647
cos3 = 0.99851488
#### aa = [0.14347029d, 3.1145815d, -0.12025561d, 0.093909436d, -0.39321884d-5, 0.10194945d-3, -0.12799464d]
#### bb = [-4.5467523d, 3.1848875d, -0.16329986d, -0.09776818d, 0.17556527d-3, -0.01978317d, 44.55915d]
aa = np.array([0.14347029, 3.1145815, -0.12025561, 0.093909436, -0.39321884e-5, 0.10194945e-3, -0.12799464])
bb = np.array([-4.5467523, 3.1848875, -0.16329986, -0.09776818, 0.17556527e-3, -0.01978317, 44.55915])
#### ; First calculate the latitude and longitude in non-rotating Jupiter coordinates
#### ; Calculate the best fit theta and fphi
#### ;t = eetime(ctime) - etime1
#### t = leapSecond(ctime) - etime1
#### cos_omeg = cos(omegay*t)
#### sin_omeg = sin(omegay*t)
t = leapSecond(ctime) - etime1
cos_omeg = np.cos(omegay * t)
sin_omeg = np.sin(omegay * t)
#### ;x = dblarr(7)
#### x5 = t/year
#### x = [cos_omeg, sin_omeg, cos(2*omegay*t), $
#### 2*cos_omeg*sin_omeg, x5*x5, x5, 1d]
x5 = t / year
x = np.array([cos_omeg, sin_omeg, np.cos(2*omegay*t), 2*cos_omeg*sin_omeg, x5*x5, x5, 1.])
#### ; fphi is phi in Jupiter fixed (non-rotating) coordinate
#### fphi = total(bb*x)
#### stheta = total(aa*x)
fphi = (bb * x).sum()
stheta = (aa * x).sum()
####
#### ; Now rotate the longitude to Jupiter System III
#### ; First Add the rotation of Jupiter around the Sun
#### ; fphi is the phi of the Sun as unspinning Jupiter goes around the sun
#### fphi = (fphi + t/yrjup*360d) mod D360
fphi = (fphi + t / yrjup * 360.) % D360
####
#### ; Next add the rotation of Jupiter around its axis.
#### sphi = (fphi - t*omega) mod D360
#### if (sphi lt 0d) then sphi = sphi + 360d
#### sphi = sphi*radian
#### stheta = stheta*radian
#### acostanstheta = acos(tan(stheta)/tan3)
sphi = (fphi - t * omega) % D360
if sphi < 0:
sphi = sphi + 360.
sphi = sphi * radian
stheta = stheta * radian
acostanstheta = np.arccos(np.tan(stheta) / tan3)
####
#### ; Now compute the orbital phase (called phi2 or phase here)
#### ; There are two solutions to the problem. Only one that is close to phi2b is correct
#### if (stheta ge three) then stheta = three
if stheta >= three:
stheta = three
#### if (-stheta ge three) then stheta = -three
if -stheta >= three:
stheta = - three
#### phi21 = (sphi + PI) + acostanstheta
phi21 = (sphi + pi) + acostanstheta
#### phi22 = (sphi + PI) - acostanstheta
phi22 = (sphi + pi) - acostanstheta
#### phi21 = phi21 mod twopi
phi21 = phi21 % twopi
#### phi22 = phi22 mod twopi
phi22 = phi22 % twopi
#### dphi = fphi + 48.23012d
dphi = fphi + 48.23012
#### phi2b = sphi - dphi*radian
phi2b = sphi - dphi * radian
#### phi2b = phi2b mod twopi
phi2b = phi2b % twopi
#### phase = phi21
phase = phi21
#### dif2 = abs(phi22 - phi2b)
dif2 = np.abs(phi22 - phi2b)
#### if (dif2 gt 350d*radian) then dif2 = twopi - dif2
if dif2 > 350. * radian:
dif2 = twopi - dif2
#### dif1 = abs(phi21 - phi2b)
dif1 = np.abs(phi21 - phi2b)
#### if (dif1 gt 350d*radian) then dif1 = twopi - dif1
if dif1 > 350. * radian:
dif1 = twopi - dif1
#### if (dif2 lt dif1) then phase = phi22
if dif2 < dif1:
phase = phi22
#### ;phase = phi2b
#### phase = phase mod twopi
phase = phase % twopi
return (stheta, sphi, phase)
#### end
####
[docs]def jrot(from_, to_, vecin, ctime):
#### ;--------------------------------------------------------------------------------------
#### ; procedure JROT
#### ;
#### ; INPUTS: From: a chracter*3 variable denoting the incoming coordinate system
#### ; To: a chracter*3 variable denoting the outgoing coordinate system
#### ; Vecin a variable of dimension 3 containing the incoming vector
#### ; Vecout: a variable of dimension 3 containing the outgoing vector
#### ; ctime: a double precision variable denoting Cline time which can be
#### ; calculated by using the function program ctimer.
#### ;
#### ; The "From" coordinate system is rotated to system III, then rotated
#### ; into the "To" coordinate system
#### ;
#### ; The supported coordinate systems are: S3C System III Cartesian (right-handed)
#### ; JSO Jupiter-Sun-Orbital
#### ; JSM Jupite-Sun-Magnetic
#### ; DIP Dipole (cartesian)
#### ; JSS Jupiter-Sun-Spin rewritten 7/2/09 by M.D.
#### ;--------------------------------------------------------------------------------------
#### pro JROT, From, To, Vecin, Vecout, ctime
''' Convert by rotation.
:param from_: Incoming frame
:param to_: Outgoing frame
:param vecin: Vector in the incoming frame.
:param ctime: a double precision variable denoting Cline time which can be
calculated by using the function program ctimer.
The "from" coordinate system is rotated to system III, then rotated into the "to" coordinate system.
The supported coordinate systems are:
* S3C System III Cartesian (right-handed)
* JSO Jupiter-Sun-Orbital
* JSM Jupite-Sun-Magnetic
* DIP Dipole (cartesian)
* JSS Jupiter-Sun-Spin
'''
#### vector = dblarr(3)
vector = np.zeros([3])
#### dipole = dblarr(3,3)
dipole = np.zeros([3, 3])
#### first = (second = dblarr(3,3) )
first = np.zeros([3, 3])
second = np.zeros([3, 3])
#### dummy = dblarr(3,3)
dummy = np.zeros([3, 3])
#### ;Identity = dblarr(3,3)
#### Identity = [ [1d, 0d, 0d], $
#### [0d, 1d, 0d], $
#### [0d, 0d, 1d] ]
identity = np.diag([1.0, 1, 1])
#### PI = 3.14159265358979d
pi = np.pi
#### twopi = 2d*PI
twopi = 2 * pi
#### radian = PI/180d
radian = pi / 180.
#### degree = 180d/PI
degree = 180. / pi
#### three = 3.123*radian
three = 3.123 * radian
#### tan3 = 0.054560676d
#### sin3 = 0.054479647d
#### cos3 = 0.99851488d
tan3 = 0.054560676
sin3 = 0.054479647
cos3 = 0.99851488
#### TH_DIP = 0.16755161d
#### PH_DIP = -3.5255651d
th_dip = 0.16755161
ph_dip = -3.5255651
#### dipole[2,*] = [-0.154625290d, 0.0624726744d, 0.98599604d]
#### dipole[1,*] = [ -0.3746066d, -0.92718385d, 0d]
#### dipole[0,*] = [ -0.9141996d, 0.36936062d, -0.16676875d]
#### ;dipole = transpose(dipole)
dipole[2,:] = [-0.154625290, 0.0624726744, 0.98599604]
dipole[1,:] = [ -0.3746066, -0.92718385, 0]
dipole[0,:] = [ -0.9141996, 0.36936062, -0.16676875]
### dipole = dipole.T #### Needed or not? TBC.
####
#### defsysv, "!stheta", exists = exist
#### if exist eq 1 then begin
#############################
# HERE, YF THINKS THE BUG. SEE WARNINGS IN THE TOP_DOCSTRING.
if _stheta != None:
#### stheta = !stheta
#### sphi = !sphi
#### phase = !phase
stheta = _stheta
sphi = _sphi
phase = _phase ## FOR ME (YF), this is extremely dangerous way of getting parameters...
## Who can guarantee the global values are as expected?
#### endif else begin
else:
#### JSun, ctime, stheta, sphi, phase
#### endelse
stheta, sphi, phase = JSun(ctime)
############################## Indeed, JSun SHOULD BE CALLED ALL THE TIME.
############################## This module keeps the bug as it is. Use kk_2009 module for nominal use.
####
#### sin_stheta = sin(stheta)
#### cos_stheta = cos(stheta)
#### sin_sphi = sin(sphi)
#### cos_sphi = cos(sphi)
sin_stheta = np.sin(stheta)
cos_stheta = np.cos(stheta)
sin_sphi = np.sin(sphi)
cos_sphi = np.cos(sphi)
####
#### case from of
#### "s3c" : first = Identity
if from_.lower() == 's3c':
first = identity
elif from_.lower() == 'jso':
#### "jso" : begin
#### ; print, " Incoming coordinate system is JSO"
#### ; Calculate the rotation matrix to go from JSO to S3R
#### ; Define X component in system III of XJSO unit vector etc.
#### dummy[0,*] = [cos_stheta*cos_sphi, cos_stheta*sin_sphi, sin_stheta]
dummy[0, :] = [cos_stheta*cos_sphi, cos_stheta*sin_sphi, sin_stheta]
#### ; Calculate the Z axis of the JSO coordinate system from the fact that the Z axis is tilted
#### ; by 3.12 degrees from the Z axis of SIII and is normal to the X axis of the JSO system.
#### ; Define X component in system III of ZJSO unit vector etc.
#### dummy[2,*] = [sin3*cos(phase), sin3*sin(phase), cos3]
dummy[2,:] = [sin3*np.cos(phase), sin3*np.sin(phase), cos3]
#### ; Define X component in system III of YJSO unit vector etc.
#### dummy[1,0] = dummy[2,1]*dummy[0,2] - dummy[2,2]*dummy[0,1]
#### dummy[1,1] = dummy[2,2]*dummy[0,0] - dummy[2,0]*dummy[0,2]
#### dummy[1,2] = dummy[2,0]*dummy[0,1] - dummy[2,1]*dummy[0,0]
dummy[1,0] = dummy[2,1]*dummy[0,2] - dummy[2,2]*dummy[0,1]
dummy[1,1] = dummy[2,2]*dummy[0,0] - dummy[2,0]*dummy[0,2]
dummy[1,2] = dummy[2,0]*dummy[0,1] - dummy[2,1]*dummy[0,0]
first = dummy.T
#### end ; jso case
elif from_.lower() == 'jsm':
#### "jsm" : begin
#### ; print, " Incoming coordinate system is JSM"
#### ; Now define the JSM transpose vector
#### ; Define X component in system III of XJSM unit vector etc.
#### dummy[0,0] = cos_stheta*cos_sphi
#### dummy[0,1] = cos_stheta*sin_sphi
#### dummy[0,2] = sin_stheta
dummy[0,0] = cos_stheta*cos_sphi
dummy[0,1] = cos_stheta*sin_sphi
dummy[0,2] = sin_stheta
#### ; Now define the Y vector so that it is perpendicular to the dipole vector and X
#### dummy[1,0] = dummy[0,2]*dipole[2,1] - dummy[0,1]*dipole[2,2]
#### dummy[1,1] = dummy[0,0]*dipole[2,2] - dummy[0,2]*dipole[2,0]
#### dummy[1,2] = dummy[0,1]*dipole[2,0] - dummy[0,0]*dipole[2,1]
dummy[1,0] = dummy[0,2]*dipole[2,1] - dummy[0,1]*dipole[2,2]
dummy[1,1] = dummy[0,0]*dipole[2,2] - dummy[0,2]*dipole[2,0]
dummy[1,2] = dummy[0,1]*dipole[2,0] - dummy[0,0]*dipole[2,1]
#### denom = sqrt(dummy[1,0]*dummy[1,0] + dummy[1,1]*dummy[1,1] + dummy[1,2]*dummy[1,2])
denom = np.sqrt(dummy[1,0]*dummy[1,0] + dummy[1,1]*dummy[1,1] + dummy[1,2]*dummy[1,2])
#### dummy[1,0] = dummy[1,0]/denom
#### dummy[1,1] = dummy[1,1]/denom
#### dummy[1,2] = dummy[1,2]/denom
dummy[1,0] = dummy[1,0]/denom
dummy[1,1] = dummy[1,1]/denom
dummy[1,2] = dummy[1,2]/denom
#### ; Now define the z vector
dummy[2,0] = dummy[0,1]*dummy[1,2] - dummy[0,2]*dummy[1,1]
dummy[2,1] = dummy[0,2]*dummy[1,0] - dummy[0,0]*dummy[1,2]
dummy[2,2] = dummy[0,0]*dummy[1,1] - dummy[0,1]*dummy[1,0]
first = dummy.T
#### end ; jsm case
elif from_.lower() == 'dip':
first = dipole.T
#### "dip" : first = transpose(dipole)
elif from_.lower() == 'jss':
#### "jss" : begin
#### ; print, " Incoming coordinate system is JSS"
#### ; Calculate the rotation matrix to go from JSS to S3R
#### ; Define X component in system III of ZJSS unit vector etc.
#### dummy[2,0] = 0d
#### dummy[2,1] = 0d
#### dummy[2,2] = 1d
dummy[2,0] = 0.
dummy[2,1] = 0.
dummy[2,2] = 1.
#### ; Define X component in system III of YJSS unit vector etc.
#### ; so that it is perpendicular to the system 3 Z spin axis and the
#### ; Jupiter-Sun line
#### dummy[1,0] = -cos_stheta*sin_sphi
#### dummy[1,1] = cos_stheta*cos_sphi
#### dummy[1,2] = 0d
dummy[1,0] = -cos_stheta*sin_sphi
dummy[1,1] = cos_stheta*cos_sphi
dummy[1,2] = 0.
denom = np.sqrt(dummy[1,0]*dummy[1,0] + dummy[1,1]*dummy[1,1])
#### dummy[1,0] = dummy[1,0]/denom
#### dummy[1,1] = dummy[1,1]/denom
#### dummy[1,2] = dummy[1,2]/denom
dummy[1,0] = dummy[1,0]/denom
dummy[1,1] = dummy[1,1]/denom
dummy[1,2] = dummy[1,2]/denom
#### ; Define X component in system III of XJSS unit vector etc.
#### ; so that it is perpendicular to ZJSS and YJSS
#### dummy[0,0] = dummy[1,1]
#### dummy[0,1] = -dummy[1,0]
#### dummy[0,2] = 0d
dummy[0,0] = dummy[1,1]
dummy[0,1] = -dummy[1,0]
dummy[0,2] = 0.
#### ; The transpose of the dummy matrix takes JSS to System 3.
#### first = transpose(dummy)
first = dummy.T
#### end ; jss case
#### endcase
####
#### case to of
if to_.lower() == 's3c':
#### "s3c" : second = identity
second = identity
elif to_.lower() == 'jso':
#### "jso" : begin
#### ; print, " The outgoing coordinate system is JSO"
#### ; Get the matrix that Rotates from System III cartesian into JSO.
#### ; Define X component in system III of XJSO unit vector etc.
#### Second[0,0] = cos_stheta*cos_sphi
#### Second[0,1] = cos_stheta*sin_sphi
#### Second[0,2] = sin_stheta
second[0,0] = cos_stheta*cos_sphi
second[0,1] = cos_stheta*sin_sphi
second[0,2] = sin_stheta
#### ; Calculate the Z axis of the JSO coordinate system from the fact that the Z axis is tilted
#### ; by 3.12 degrees from the Z axis of SIII and is normal to the X axis of the JSO system.
#### ; Define X component in system III of ZJSO unit vector etc.
#### Second[2,0] = sin3*cos(phase)
#### Second[2,1] = sin3*sin(phase)
#### Second[2,2] = cos3
second[2,0] = sin3*np.cos(phase)
second[2,1] = sin3*np.sin(phase)
second[2,2] = cos3
#### ; Define X component in system III of YJSO unit vector etc.
#### Second[1,0] = Second[2,1]*Second[0,2] - Second[2,2]*Second[0,1]
#### Second[1,1] = Second[2,2]*Second[0,0] - Second[2,0]*Second[0,2]
#### Second[1,2] = Second[2,0]*Second[0,1] - Second[2,1]*Second[0,0]
second[1,0] = second[2,1]*second[0,2] - second[2,2]*second[0,1]
second[1,1] = second[2,2]*second[0,0] - second[2,0]*second[0,2]
second[1,2] = second[2,0]*second[0,1] - second[2,1]*second[0,0]
#### end ; jso case
elif to_.lower() == 'jsm':
#### "jsm" : begin
#### ; print, " The outgoing coordinate system is JSM"
#### ; Get the matrix that Rotates from System III cartesian into JSM.
#### ; Define X component in system III of XJSM unit vector etc.
#### Second[0,0] = cos_stheta*cos_sphi
#### Second[0,1] = cos_stheta*sin_sphi
#### Second[0,2] = sin_stheta
second[0,0] = cos_stheta*cos_sphi
second[0,1] = cos_stheta*sin_sphi
second[0,2] = sin_stheta
#### ; Now define the Y vector so that it is perpendicular to the dipole vector and X
#### Second[1,0] = second[0,2]*dipole[2,1] - second[0,1]*dipole[2,2]
#### Second[1,1] = second[0,0]*dipole[2,2] - second[0,2]*dipole[2,0]
#### Second[1,2] = second[0,1]*dipole[2,0] - second[0,0]*dipole[2,1]
second[1,0] = second[0,2]*dipole[2,1] - second[0,1]*dipole[2,2]
second[1,1] = second[0,0]*dipole[2,2] - second[0,2]*dipole[2,0]
second[1,2] = second[0,1]*dipole[2,0] - second[0,0]*dipole[2,1]
#### denom = sqrt(second[1,0]*second[1,0] + second[1,1]*second[1,1] + second[1,2]*second[1,2])
denom = np.sqrt(second[1,0]*second[1,0] + second[1,1]*second[1,1] + second[1,2]*second[1,2])
#### Second[1,0] = second[1,0]/denom
#### Second[1,1] = second[1,1]/denom
#### Second[1,2] = second[1,2]/denom
second[1,0] = second[1,0]/denom
second[1,1] = second[1,1]/denom
second[1,2] = second[1,2]/denom
#### ; Now define the z vector
#### Second[2,0] = second[0,1]*second[1,2] - second[0,2]*second[1,1]
#### Second[2,1] = second[0,2]*second[1,0] - second[0,0]*second[1,2]
#### Second[2,2] = second[0,0]*second[1,1] - second[0,1]*second[1,0]
second[2,0] = second[0,1]*second[1,2] - second[0,2]*second[1,1]
second[2,1] = second[0,2]*second[1,0] - second[0,0]*second[1,2]
second[2,2] = second[0,0]*second[1,1] - second[0,1]*second[1,0]
#### end ; jsm case
elif to_.lower() == 'dip':
second = dipole.T ### YF considers this may be a bug, or the other one (first)...
#### "dip" : second = transpose(dipole)
elif to_.lower() == 'jss':
#### "jss" : begin
#### ; print, " The outgoing coordinate system is JSS"
#### ; Define X component in system III of ZJSS unit vector etc.
#### second[2,0] = 0d
#### second[2,1] = 0d
#### second[2,2] = 1d
second[2,0] = 0.
second[2,1] = 0.
second[2,2] = 1.
#### ; Define X component in system III of YJSS unit vector etc.
#### ; so that it is perpendicular to the system 3 Z spin axis and the
#### ; Jupiter-Sun line
#### second[1,0] = -cos_stheta*sin_sphi
#### second[1,1] = cos_stheta*cos_sphi
#### second[1,2] = 0d
second[1,0] = -cos_stheta*sin_sphi
second[1,1] = cos_stheta*cos_sphi
second[1,2] = 0.
#### denom = sqrt(second[1,0]*second[1,0] + second[1,1]*second[1,1] + second[1,2]*second[1,2])
denom = np.sqrt(second[1,0]*second[1,0] + second[1,1]*second[1,1] + second[1,2]*second[1,2])
#### second[1,0] = second[1,0]/denom
#### second[1,1] = second[1,1]/denom
#### second[1,2] = second[1,2]/denom
second[1,0] = second[1,0]/denom
second[1,1] = second[1,1]/denom
second[1,2] = second[1,2]/denom
#### ; Define X component in system III of XJSS unit vector etc.
#### ; so that it is perpendicular to ZJSS and YJSS
#### second[0,0] = second[1,1]*second[2,2] - second[1,2]*second[2,1]
#### second[0,1] = second[1,2]*second[2,0] - second[1,0]*second[2,2]
#### second[0,2] = second[1,0]*second[2,1] - second[1,1]*second[2,0]
second[0,0] = second[1,1]*second[2,2] - second[1,2]*second[2,1]
second[0,1] = second[1,2]*second[2,0] - second[1,0]*second[2,2]
second[0,2] = second[1,0]*second[2,1] - second[1,1]*second[2,0]
#### end ; jss case
#### endcase
####
#### ; Now multimply vecin with first and second matrices to get the vecout
#### vector = matrix_multiply(first, vecin)
vector = first.dot(vecin)
#### vecout = matrix_multiply(second, vector)
vecout = second.dot(vector)
return vecout
#### end ; jrot procedure
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure cy12car_pos
#### ;--------------------------------------------------------------------------------------
#### pro cyl2car_pos, rho, phi, Z, Xcar, Ycar, Zcar
[docs]def cyl2car_pos(rho, phi, Z):
#### Xcar = rho*cos(phi)
#### Ycar = rho*sin(phi)
#### Zcar = Z
Xcar = rho*np.cos(phi)
Ycar = rho*np.sin(phi)
Zcar = Z
return (Xcar, Ycar, Zcar)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure car2cyl_pos
#### ;--------------------------------------------------------------------------------------
#### pro car2cyl_pos, rho, phi, Z, Xcar, Ycar, Zcar
#### phi = atan(Ycar,Xcar)
#### rho = sqrt(Xcar*Xcar + Ycar*Ycar)
#### Z = Zcar
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure car2cyl_vect
#### ;--------------------------------------------------------------------------------------
#### pro car2cyl_vect, Brho, Bphi, BZ, BXcar, BYcar, BZcar, phi
[docs]def car2cyl_vect(BXcar, BYcar, BZcar, phi):
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
Brho = BXcar*cos_phi + BYcar*sin_phi
Bphi = -BXcar*sin_phi + BYcar*cos_phi
BZ = BZcar
return (Brho, Bphi, BZ)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure cy12car_vect
#### ;--------------------------------------------------------------------------------------
#### pro cyl2car_vect, Brho, Bphi, BZ, BXcar, BYcar, BZcar, phi
[docs]def cyl2car_vect(Brho, Bphi, BZ, phi):
#### cos_phi = cos(phi)
#### sin_phi = sin(phi)
#### BXcar = Brho*cos_phi - Bphi*sin_phi
#### BYcar = Brho*sin_phi + Bphi*cos_phi
#### BZcar = BZ
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
BXcar = Brho*cos_phi - Bphi*sin_phi
BYcar = Brho*sin_phi + Bphi*cos_phi
BZcar = BZ
return BXcar, BYcar, BZcar
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure CAR2SPH_pos
#### ;--------------------------------------------------------------------------------------
#### pro CAR2SPH_pos, X, Y, Z, R, TH, PHI
[docs]def CAR2SPH_pos(X, Y, Z):
#### ; BARTSCH
#### R = sqrt(X*X + Y*Y + Z*Z)
#### PHI = atan(Y, X)
#### TH = asin(Z/R)
R = np.sqrt(X*X + Y*Y + Z*Z)
PHI = np.arctan2(Y, X)
TH = np.arcsin(Z/R)
return (R, TH, PHI)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure SPH2CAR_pos
#### ;--------------------------------------------------------------------------------------
#### pro SPH2CAR_pos, X, Y, Z, R, TH, PHI
[docs]def SPH2CAR_pos(R, TH, PHI):
#### ; BARTSCH
#### sin_th = sin(TH)
#### X = R*sin_th*cos(PHI)
#### Y = R*sin_th*sin(PHI)
#### Z = R*cos(TH)
sin_th = np.sin(TH)
X = R*sin_th*np.cos(PHI)
Y = R*sin_th*np.sin(PHI)
Z = R*np.cos(TH)
return (X, Y, Z)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure SPH2CAR_MAG
#### ;--------------------------------------------------------------------------------------
#### pro SPH2CAR_MAG, BX, BY, BZ, BR, BTH, BPHI, TH, PHI
#### ; ARFKEN
#### sin_th = sin(TH)
#### cos_th = cos(TH)
#### sin_phi = sin(PHI)
#### cos_phi = cos(PHI)
#### BX = BR*sin_th*cos_phi + BTH*cos_th*cos_phi - BPHI*sin_phi
#### BY = BR*sin_th*sin_phi + BTH*cos_th*sin_phi + BPHI*cos_phi
#### BZ = BR*cos_th - BTH*sin_th
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure CAR2SPH_MAG
#### ;--------------------------------------------------------------------------------------
#### pro CAR2SPH_MAG, BX, BY, BZ, BR, BTH, BPHI, TH, PHI
[docs]def CAR2SPH_MAG(BX, BY, BZ, TH, PHI):
#### ; ARFKEN
#### sin_th = sin(TH)
sin_th = np.sin(TH)
#### cos_th = cos(TH)
cos_th = np.cos(TH)
#### sin_phi = sin(PHI)
sin_phi = np.sin(PHI)
#### cos_phi = cos(PHI)
cos_phi = np.cos(PHI)
#### BR = BX*sin_th*cos_phi + BY*sin_th*sin_phi + BZ*cos_th
BR = BX*sin_th*cos_phi + BY*sin_th*sin_phi + BZ*cos_th
#### BTH = BX*cos_th*cos_phi + BY*cos_th*sin_phi - BZ*sin_th
BTH = BX*cos_th*cos_phi + BY*cos_th*sin_phi - BZ*sin_th
#### BPHI = -BX*sin_phi + BY*cos_phi
BPHI = -BX*sin_phi + BY*cos_phi
return (BR, BTH, BPHI)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure PJSO2S3
#### ;--------------------------------------------------------------------------------------
#### pro PJSO2S3, BXPJSO, BYPJSO, BZPJSO, BXS3, BYS3, BZS3, phi
[docs]def PJSO2S3(BXPJSO, BYPJSO, BZPJSO, phi):
sin_phi = np.sin(phi)
cos_phi = np.cos(phi)
BXS3 = BXPJSO*cos_phi - BYPJSO*sin_phi
BYS3 = BXPJSO*sin_phi + BYPJSO*cos_phi
BZS3 = BZPJSO
return (BXS3, BYS3, BZS3)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure S32PJSO
#### ;--------------------------------------------------------------------------------------
#### pro S32PJSO, BXPJSO, BYPJSO, BZPJSO, BXS3, BYS3, BZS3, phi
[docs]def S32PJSO(BXS3, BYS3, BZS3, phi):
#### sin_phi = sin(phi)
#### cos_phi = cos(phi)
#### BXPJSO = BXS3*cos_phi + BYS3*sin_phi
#### BYPJSO = -BXS3*sin_phi + BYS3*cos_phi
#### BZPJSO = BZS3
sin_phi = np.sin(phi)
cos_phi = np.cos(phi)
BXPJSO = BXS3*cos_phi + BYS3*sin_phi
BYPJSO = -BXS3*sin_phi + BYS3*cos_phi
BZPJSO = BZS3
return (BXPJSO, BYPJSO, BZPJSO)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; function bessj1(x)
#### ; Taken from: http://audiolab.uwaterloo.ca/~jeffb/thesis/node50.html
#### ;--------------------------------------------------------------------------------------
[docs]def bessj1(x):
''' Bessel function.
Returns bessel function with real order of 1.
>>> print('{:.6f}'.format(bessj1(0)))
0.000000
>>> print('{:.6f}'.format(bessj1(0.1)))
0.049938
>>> print('{:.6f}'.format(bessj1(2)))
0.576725
>>> print('{:.6f}'.format(bessj1(-0.1)))
-0.049938
Internally, :func:`scipy:scipy.special.jv` is used.
The returned values agreed with the original IDL version.
'''
return special.jv(1, x)
### function bessj1, x
#### a = abs([min(x), max(x)])
#### if (a[0] gt a[1]) then b = [a[1],a[0]] else b = a
#### ; everything here is between -8 and 8
#### if (b[1] lt 8d) then begin
#### y = x*x
#### bessj1 = x*(72362614232d + y*(-7895059235d + y*(242396853.1d + y*(-2972611.439d + y*( 15704.48260d + y*(-30.160366606d) ))))) $
#### /(144725228442d + y*(2300535178d + y*(18583304.74d + y*(99447.43394d + y*(376.9991397d + y*1.0d)))))
#### return, bessj1
#### endif
#### ; everything is either < or = to -8 or > or = to +8
#### if (b[0] ge 8d) then begin
#### ax = abs(x)
#### z = 8d/ax
#### y = z*z
#### xx = ax - 2.356194491d
#### if x ge 0d then one = 1d else one = -1d
#### bessj1 = sqrt(0.636619772d/ax)*(cos(xx)*(1.0d + y*(0.183105d-2 + y*(-0.3516396496d-4 + y*( 0.2457520174d-5 + y*(-0.240337019d-6) )))) $
#### -z*sin(xx)*(0.04687499995d0 + y*(-0.88228987d-6 + y*(0.105787412d-6 + y*(-0.88228987d-6 + y*0.105787412d-6)))))*one
#### return, bessj1
#### endif
#### ; below is a mixture of the two cases above
#### y = x
#### for ii = 0, n_elements(x) - 1 do y[ii] = bessj1(x[ii])
#### bessj1 = y
#### return, bessj1
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; function bessj0(x)
#### ; Taken from: http://audiolab.uwaterloo.ca/~jeffb/thesis/node50.html
#### ;--------------------------------------------------------------------------------------
[docs]def bessj0(x):
''' bessel function with real order of 0.
>>> print('{:.6f}'.format(bessj0(0)))
1.000000
>>> print('{:.6f}'.format(bessj0(0.1)))
0.997502
>>> print('{:.6f}'.format(bessj0(2)))
0.223891
>>> print('{:.6f}'.format(bessj0(-0.1)))
0.997502
Internally, :func:`scipy:scipy.special.jv` is used.
The returned values agreed with the original IDL version.
'''
return special.jv(0, x)
#### function bessj0, x
#### a = abs([min(x), max(x)])
#### if (a[0] gt a[1]) then b = [a[1],a[0]] else b = a
#### ; everything here is between -8 and 8
#### if (b[1] lt 8) then begin
#### y = x*x
#### bessj0 = (57568490574d + y*(-13362590354d + y*(651619640.7d + y*(-11214424.18d + y*( 77392.33017d - 184.9052456d*y))))) $
#### /(57568490411d + y*(1029532985d + y*(9494680.718d + y*(59272.64853d + y*(267.8532712d + y*1.0d)))))
#### return, bessj0
#### endif
#### ; everything is either < or = to -8 or > or = to +8
#### if (b[0] ge 8) then begin
#### ax = abs(x)
#### z = 8d/ax
#### y = z*z
#### xx = ax - 0.785398164d
#### bessj0 = sqrt(0.636619772d/ax)*(cos(xx)*(1.0d + y*(-.1098628627d-2 + y*(.2734510407d-4 + y*(-.2073370639d-5 + y*.2093887211d-6)))) $
#### - z*sin(xx)*(-.1562499995d-1 + y*(.1430488765d-3 + y*(-.6911147651d-5 + y*(.7621095161d-6 - .934945152d-7*y)))))
#### return, bessj0
#### endif
#### ; below is a mixture of the two cases above
#### y = x
#### for ii = 0, n_elements(x) - 1 do y[ii] = bessj0(x[ii])
#### bessj0 = y
#### return, bessj0
####
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; function DBSJ2(x)
#### ;--------------------------------------------------------------------------------------
#### function DBSJ2, x
[docs]def DBSJ2(x):
return (2. / x) * bessj1(x) - bessj0(x)
#### DBSJ2 = (2d/x)*bessj1(x) - bessj0(x)
#### return, DBSJ2
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; function DBSJ3(x)
#### ;--------------------------------------------------------------------------------------
#### function DBSJ3, x
[docs]def DBSJ3(x):
#### DBSJ3 = (4d/x)*DBSJ2(x) - bessj1(x)
return (4./x)*DBSJ2(x) - bessj1(x)
#### ;print, (4d/x)*DBSJ2(x) - bessj1(x)
#### ;print, (8d/(x*x) - 1)*bessj1(x) - (4d/x)*bessj0(x)
#### ;print, ''
#### ;DBSJ3 = (8d/(x*x) - 1)*bessj1(x) - (4d/x)*bessj0(x) ; RJW
#### return, DBSJ3
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure B_mp_perp
#### ;--------------------------------------------------------------------------------------
#### pro B_mp_perp, rho, phi, x, Bperpr, Bperpf, Bperpx, Nmodes
[docs]def B_mp_perp(rho, phi, x, Nmodes):
#### IDL sample.pro (first call) returns
# rho = 5.7614139
# phi = 1.5707963
# x = 29.774504
# Nmodes = 3
#### Python sreturns
# >>> rho, phi, x, Nmodes
# (5.7614152250836614, 1.5707963267948966, 29.774497437565074, 3)
# So far so good 2014-05-04T12:32:35
#### ; A
#### A = [14.8636286417325d, -20.6733061430101, -68.9799158183312d, -1.42668213471797D-006, $
#### -2.10151726837109D-004, 7.79886406933501D-006, 1.68914984127525D-008, -1.07213724493682D-007, $
#### 5.22885344671751D-007, -4.46205598371400D-008, 2.83331673403534D-007, -1.38338264037040D-006, $
#### 7.94736850421990d, 23.4906454129448d, -17.3725909963207d, -3.95338621079697D-004, $
#### -4.24194742905271D-003, 1.92089873646750D-003, -6.94191156816239D-002, -0.170169494909972d, $
#### -0.171632981837208d, 4.57144670517197D-009, 6.53500870745233D-007, 6.56440399411187D-007, $
#### -7.67358819268301D-002, -0.133394608915796d, -0.127893292811889d, 3.60644094744725D-006 , $
#### 2.04677410835679D-005, 3.26595444981956D-005]
####
A = np.array([14.8636286417325, -20.6733061430101, -68.9799158183312, -1.42668213471797e-6,
-2.10151726837109e-4, 7.79886406933501e-6, 1.68914984127525e-8, -1.07213724493682e-7,
5.22885344671751e-7, -4.46205598371400e-8, 2.83331673403534e-7, -1.38338264037040e-6,
7.94736850421990, 23.4906454129448, -17.3725909963207, -3.95338621079697e-4,
-4.24194742905271e-3, 1.92089873646750e-3, -6.94191156816239e-2, -0.170169494909972,
-0.171632981837208, 4.57144670517197e-9, 6.53500870745233e-7, 6.56440399411187e-7,
-7.67358819268301e-2, -0.133394608915796, -0.127893292811889, 3.60644094744725e-6 ,
2.04677410835679e-5, 3.26595444981956e-5])
#### ; the b(J) are
#### B = [32.9198722839323d, 64.1556701660173d, 128.854949951183d, 30.4416656494141d, $
#### 45.2417564392090d, 107.303962707520d, 28.7235050200276d, 63.9741210937499d, $
#### 128.002441406250d, 33.8896331786845d, 60.2880744933923d, 119.902839660605d, $
#### 32.1036300655818d, 62.8152580261215d, 127.073867797851d]
B = np.array([32.9198722839323, 64.1556701660173, 128.854949951183, 30.4416656494141,
45.2417564392090, 107.303962707520, 28.7235050200276, 63.9741210937499,
128.002441406250, 33.8896331786845, 60.2880744933923, 119.902839660605,
32.1036300655818, 62.8152580261215, 127.073867797851])
####
#### Modes1 = Nmodes
#### Modes2 = 2*Modes1
#### Modes3 = 3*Modes1
#### Modes4 = 4*Modes1
#### Modes5 = 5*Modes1
#### ;Modes6 = 6*Modes1 ; not used
####
Modes1 = Nmodes
Modes2 = 2*Modes1
Modes3 = 3*Modes1
Modes4 = 4*Modes1
Modes5 = 5*Modes1
####
#### rhom = rho - 0.1d
#### rhop = rho + 0.1d
####
rhom = rho - 0.1
rhop = rho + 0.1
#### ;phim = phi - 0.001d
#### ;phip = phi + 0.001d
####
#### phi2 = 2d*phi
phi2 = 2.*phi
#### ;phi2m = 2d*phi - 0.001d
#### ;phi2p = 2d*phi + 0.001d
####
#### phi3 = 3d*phi
phi3 = 3.*phi
#### ;phi3m = 3d*phi - 0.001d
#### ;phi3p = 3d*phi + 0.001d
####
#### xm = x - 0.1d
#### xp = x + 0.1d
xm = x - 0.1
xp = x + 0.1
####
#### sinf = sin(phi)
#### cosf = cos(phi)
sinf = np.sin(phi)
cosf = np.cos(phi)
#### ;sin2f = sin(phi2)
#### sin2f = 2*sinf*cosf ; double angle
sin2f = 2*sinf*cosf
#### ;cos2f = cos(phi2)
#### cos2f = 2*cosf*cosf - 1 ; double angle
cos2f = 2*cosf*cosf - 1
#### sin3f = sin(phi3)
#### cos3f = cos(phi3)
sin3f = np.sin(phi3)
cos3f = np.cos(phi3)
####
#### expb = exp(x/b)
expb = np.exp(x/B)
#### exp_xm_minus_xp = exp(xp/b) - exp(xm/b)
exp_xm_minus_xp = np.exp(xp/B) - np.exp(xm/B)
####
#### sin_001 = sin(0.001d)
sin_001 = np.sin(0.001)
####
#### ;sinP_minus_sinM = sin(phip) - sin(phim)
#### sinP_minus_sinM = 2d*cos(phi)*sin_001
sinP_minus_sinM = 2.*np.cos(phi)*sin_001
####
#### ;cosP_minus_cosM = cos(phip) - cos(phim)
#### cosP_minus_cosM = -2d*sin(phi)*sin_001
cosP_minus_cosM = -2.*np.sin(phi)*sin_001
####
#### ;sin2P_minus_sin2M = sin(phi2p) - sin(phi2m)
#### sin2P_minus_sin2M = 2d*cos(phi2)*sin_001
sin2P_minus_sin2M = 2.*np.cos(phi2)*sin_001
####
#### ;cos2P_minus_cos2M = cos(phi2p) - cos(phi2m)
#### cos2P_minus_cos2M = -2d*sin(phi2)*sin_001
cos2P_minus_cos2M = -2.*np.sin(phi2)*sin_001
####
#### ;sin3P_minus_sin3M = sin(phi3p) - sin(phi3m)
#### sin3P_minus_sin3M = 2d*cos(phi3)*sin_001
sin3P_minus_sin3M = 2.*np.cos(phi3)*sin_001
####
#### ;cos3P_minus_cos3M = cos(phi3p) - cos(phi3m)
#### cos3P_minus_cos3M = -2d*sin(phi3)*sin_001
cos3P_minus_cos3M = -2.*np.sin(phi3)*sin_001
####
#### rhobb = rho/b
rhobb = rho/B
#### xmb = (x - b)
xmb = (x - B)
#### fiveHundredoverRho = (500d/rho)
fiveHundredoverRho = (500./rho)
####
#### Bperpr = 0d
#### Bperpf = 0d
#### Bperpx = 0d
#### KA = 0
Bperpr = 0.
Bperpf = 0.
Bperpx = 0.
KA = 0
####
#### ;--------------------------------------------------------------------------------------
#### ; Now compute the terms associated with phi
#### Nmodes_indgen = indgen(Nmodes)
#### Nmodes_x2 = 2*Nmodes
Nmodes_indgen = np.arange(Nmodes)
Nmodes_x2 = 2*Nmodes
####
#### ind = [Nmodes_indgen]
ind = Nmodes_indgen
#### KA1 = KA + Nmodes_indgen ; KA = 0 on entering here - just here for clarity with other bits below
KA1 = KA + Nmodes_indgen
#### KA_Temp = ind + Nmodes
KA_Temp = ind + Nmodes
#### bessj1PminusM = bessj1(rhop/b[ind]) - bessj1(rhom/b[ind])
bessj1PminusM = bessj1(rhop/B[ind]) - bessj1(rhom/B[ind])
#### bessj1RHO = bessj1(rhobb[ind])
bessj1RHO = bessj1(rhobb[ind])
#### angle_diff = sinf*a[KA1] + cosf*a[KA_Temp]
angle_diff = sinf*A[KA1] + cosf*A[KA_Temp]
####
#### Bperpr += expb[ind]*bessj1PminusM *angle_diff
#### Bperpf += expb[ind]*bessj1RHO *(sinP_minus_sinM*a[KA1] + cosP_minus_cosM*a[KA_Temp])
#### Bperpx += bessj1RHO*exp_xm_minus_xp[ind]*angle_diff
#### KA += 2*Nmodes
Bperpr += expb[ind]*bessj1PminusM *angle_diff
Bperpf += expb[ind]*bessj1RHO *(sinP_minus_sinM*A[KA1] + cosP_minus_cosM*A[KA_Temp])
Bperpx += bessj1RHO*exp_xm_minus_xp[ind]*angle_diff
KA += 2*Nmodes
####
#### ; Now compute the terms associated with 2*phi
#### ind = [Nmodes_indgen + Nmodes]
ind = Nmodes_indgen + Nmodes
#### KA1 = KA + Nmodes_indgen
KA1 = KA + Nmodes_indgen
#### KA_Temp = KA1 + Nmodes
KA_Temp = KA1 + Nmodes
#### dbsj2PminusM = DBSJ2(rhop/b[ind]) - DBSJ2(rhom/b[ind])
#### dbsj2Rhobb = DBSJ2(rhobb[ind])
#### angle_diff = sin2f*a[KA1] + cos2f*a[KA_Temp]
dbsj2PminusM = DBSJ2(rhop/B[ind]) - DBSJ2(rhom/B[ind])
dbsj2Rhobb = DBSJ2(rhobb[ind])
angle_diff = sin2f*A[KA1] + cos2f*A[KA_Temp]
####
#### Bperpr += expb[ind] *dbsj2PminusM *angle_diff
#### Bperpf += expb[ind] *dbsj2Rhobb *(sin2P_minus_sin2M*a[KA1] + cos2P_minus_cos2M*a[KA_Temp])
#### Bperpx += dbsj2Rhobb*exp_xm_minus_xp[ind]*angle_diff
#### KA += Nmodes_x2
Bperpr += expb[ind] *dbsj2PminusM *angle_diff
Bperpf += expb[ind] *dbsj2Rhobb *(sin2P_minus_sin2M*A[KA1] + cos2P_minus_cos2M*A[KA_Temp])
Bperpx += dbsj2Rhobb*exp_xm_minus_xp[ind]*angle_diff
KA += Nmodes_x2
####
#### ; Now compute the terms associated with 3*phi
#### ind = [Nmodes_indgen + Modes2]
ind = Nmodes_indgen + Modes2
#### KA1 = KA + Nmodes_indgen
KA1 = KA + Nmodes_indgen
#### KA_Temp = KA1 + Nmodes
KA_Temp = KA1 + Nmodes
#### dbsj3PminusM = DBSJ3(rhop/b[ind]) - DBSJ3(rhom/b[ind])
#### dbsj3Rhobb = DBSJ3(rhobb[ind])
#### angle_diff = sin3f*a[KA1] + cos3f*a[KA_Temp]
dbsj3PminusM = DBSJ3(rhop/B[ind]) - DBSJ3(rhom/B[ind])
dbsj3Rhobb = DBSJ3(rhobb[ind])
angle_diff = sin3f*A[KA1] + cos3f*A[KA_Temp]
####
#### Bperpr += expb[ind] *dbsj3PminusM *angle_diff
#### Bperpf += expb[ind] *dbsj3Rhobb *(sin3P_minus_sin3M*a[KA1] + cos3P_minus_cos3M*a[KA_Temp])
#### Bperpx += dbsj3Rhobb*exp_xm_minus_xp[ind]*angle_diff
#### KA += Nmodes_x2
Bperpr += expb[ind] *dbsj3PminusM *angle_diff
Bperpf += expb[ind] *dbsj3Rhobb *(sin3P_minus_sin3M*A[KA1] + cos3P_minus_cos3M*A[KA_Temp])
Bperpx += dbsj3Rhobb*exp_xm_minus_xp[ind]*angle_diff
KA += Nmodes_x2
####
#### ; Now include the terms associated with the derivative term
#### ind = [Nmodes_indgen + Modes3]
ind = Nmodes_indgen + Modes3
#### KA1 = KA + Nmodes_indgen
KA1 = KA + Nmodes_indgen
#### KA_Temp = KA1 + Nmodes
KA_Temp = KA1 + Nmodes
#### besj0 = bessj0(rhobb[ind])
#### besj1 = bessj1(rhobb[ind])
besj0 = bessj0(rhobb[ind])
besj1 = bessj1(rhobb[ind])
#### BJ0minusBJ1 = rhop*bessj0(rhop/b[ind]) - rhom*bessj0(rhom/b[ind]) - xmb[ind]*(bessj1(rhom/b[ind]) - bessj1(rhop/b[ind]))
BJ0minusBJ1 = rhop*bessj0(rhop/B[ind]) - rhom*bessj0(rhom/B[ind]) - xmb[ind]*(bessj1(rhom/B[ind]) - bessj1(rhop/B[ind]))
#### expxp_minus_expxm = exp(xp/b[ind])*(rho*besj0 + (xp - b[ind])*besj1) - exp(xm/b[ind])*(rho*besj0 + (xm - b[ind])*besj1)
expxp_minus_expxm = np.exp(xp/B[ind])*(rho*besj0 + (xp - B[ind])*besj1) - np.exp(xm/B[ind])*(rho*besj0 + (xm - B[ind])*besj1)
#### BJ0plusBJ1 = rho*besj0 + xmb[ind]*besj1
BJ0plusBJ1 = rho*besj0 + xmb[ind]*besj1
#### angle_diff = sinf*a[KA1] + cosf*a[KA_Temp]
angle_diff = sinf*A[KA1] + cosf*A[KA_Temp]
####
#### Bperpr += expb[ind] *BJ0minusBJ1 *angle_diff
Bperpr += expb[ind] *BJ0minusBJ1 *angle_diff
#### Bperpf += expb[ind] *BJ0plusBJ1 *(sinP_minus_sinM*a[KA1] + cosP_minus_cosM*a[KA_Temp])
Bperpf += expb[ind] *BJ0plusBJ1 *(sinP_minus_sinM*A[KA1] + cosP_minus_cosM*A[KA_Temp])
#### Bperpx += expxp_minus_expxm* angle_diff
Bperpx += expxp_minus_expxm* angle_diff
#### KA += Nmodes_x2
KA += Nmodes_x2
####
#### ind = [Nmodes_indgen + Modes4]
ind = Nmodes_indgen + Modes4
#### KA1 = KA + Nmodes_indgen
KA1 = KA + Nmodes_indgen
#### KA_Temp = KA1 + Nmodes
KA_Temp = KA1 + Nmodes
#### besj2 = DBSJ2(rhobb[ind])
besj2 = DBSJ2(rhobb[ind])
#### besj3 = DBSJ3(rhobb[ind])
besj3 = DBSJ3(rhobb[ind])
#### dbsj2minusdbsj3 = rhop*DBSJ2(rhop/b[ind]) - rhom*DBSJ2(rhom/b[ind]) - (x - 3*b[ind])*(DBSJ3(rhom/b[ind]) - DBSJ3(rhop/b[ind]))
dbsj2minusdbsj3 = rhop*DBSJ2(rhop/B[ind]) - rhom*DBSJ2(rhom/B[ind]) - (x - 3*B[ind])*(DBSJ3(rhom/B[ind]) - DBSJ3(rhop/B[ind]))
#### expxp_minus_expxm = exp(xp/b[ind])*(rho*besj2 + (xp - 3*b[ind])*besj3) - exp(xm/b[ind])*(rho*besj2 + (xm - 3*b[ind])*besj3)
expxp_minus_expxm = np.exp(xp/B[ind])*(rho*besj2 + (xp - 3*B[ind])*besj3) - np.exp(xm/B[ind])*(rho*besj2 + (xm - 3*B[ind])*besj3)
#### BJ2minusBJ3 = rho*besj2 + (x - 3*b[ind])*besj3
BJ2minusBJ3 = rho*besj2 + (x - 3*B[ind])*besj3
#### angle_diff = sin3f*a[KA1] + cos3f*a[KA_Temp]
angle_diff = sin3f*A[KA1] + cos3f*A[KA_Temp]
####
#### Bperpr += expb[ind] *dbsj2minusdbsj3*angle_diff
#### Bperpf += expb[ind] *BJ2minusBJ3 *(sin3P_minus_sin3M*a[KA1] + cos3P_minus_cos3M*a[KA_Temp])
#### Bperpx += expxp_minus_expxm *angle_diff
Bperpr += expb[ind] *dbsj2minusdbsj3*angle_diff
Bperpf += expb[ind] *BJ2minusBJ3 *(sin3P_minus_sin3M*A[KA1] + cos3P_minus_cos3M*A[KA_Temp])
Bperpx += expxp_minus_expxm *angle_diff
#### ;--------------------------------------------------------------------------------------
####
#### Bperpr = total(Bperpr) * 5d
#### Bperpf = total(Bperpf) * fiveHundredoverRho
#### Bperpx = total(Bperpx) * 5d
Bperpr = Bperpr.sum() * 5.
Bperpf = Bperpf.sum() * fiveHundredoverRho
Bperpx = Bperpx.sum() * 5.
#### IDL sample.pro (first call) returns
# Bperpr, Bperpf, Bperpx
# -0.53848543 1.2894769e-06 -0.026515011
#### Python sreturns
# -0.5384820597983946, 1.2850552969994733e-06, -0.026515192286147078
### So far so good.
_logger.debug('Bperpr, Bperpf, Bperpx = {}, {}, {}'.format(
Bperpr, Bperpf, Bperpx))
return (Bperpr, Bperpf, Bperpx)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure dipole_cyl
#### ; not used
#### ;--------------------------------------------------------------------------------------
#### pro dipole_cyl, rmap, pmap, zmap, Brd, Bpd, Bzd, counter, ctime
#### vecin = (vecou = dblarr(3))
####
#### cyl2car_pos, rmap, pmap, zmap, Xmap, Ymap, Zmap
####
#### ; convert to jso
#### vecin = [Xmap, Ymap, Zmap]
#### JROT, "s3c", "jso", vecin, vecou, ctime
#### B0x = 0d
#### B0y = 0d
#### zz = 420543d*420543d + 65920d*65920d + 24992d*24992d
#### B0z = sqrt(zz) ;Vip4 dipole magnitude
####
#### dipole, B0x, B0y, B0z, vecou[0], vecou[1], vecou[2], Bx, By, Bz
####
#### ; convert back to s3c
#### vecin = [BX, BY, BZ]
#### JROT, "jso", "s3c", vecin, vecou, ctime
#### car2cyl_vect, Brd, Bpd, Bzd, vecou[0], vecou[1], vecou[2], pmap
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure dipole
#### ;
#### ; Calculates the field of Jupiter's dipole for shielding in the magnetopause
#### ;
#### ; (B0x, B0y, B0z) is the dipole moment
#### ; x, y, z is the position vector
#### ; Bx, By, Bz is the output field vector
#### ;--------------------------------------------------------------------------------------
#### pro dipole, B0x, B0y, B0z, x, y ,z, Bx, By, Bz
[docs]def dipole(B0x, B0y, B0z, x, y, z):
''' Dipole
::
Original comments in IDL
; Calculates the field of Jupiter's dipole for shielding in the magnetopause
;
; (B0x, B0y, B0z) is the dipole moment
; x, y, z is the position vector
; Bx, By, Bz is the output field vector
'''
#### r = sqrt(x*x + y*y + z*z)
r = np.sqrt(x*x + y*y + z*z)
#### ;a = dblarr(3,3)
#### a = [ [3*x*x - r*r, 3*x*y, 3*x*z], $
#### [3*x*y, 3*y*y - r*r, 3*y*z], $
#### [3*x*z, 3*y*z, 3*z*z - r*r] ]
a = np.array([ [3*x*x - r*r, 3*x*y, 3*x*z],
[3*x*y, 3*y*y - r*r, 3*y*z],
[3*x*z, 3*y*z, 3*z*z - r*r] ])
#### r5 = r*r*r*r*r
r5 = r*r*r*r*r
#### a = a/r5
a = a/r5
Bx = a[0,0]*B0x + a[0,1]*B0y + a[0,2]*B0z
By = a[1,0]*B0x + a[1,1]*B0y + a[1,2]*B0z
Bz = a[2,0]*B0x + a[2,1]*B0y + a[2,2]*B0z
return (Bx, By, Bz)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure dipole_shielded
#### ;
#### ; WRITTEN BY K. K. KHURANA 11/2002
#### ; MODIFIED BY H. K. SCHWARZL 11/2003
#### ; PARMOD is an input array (real*8) that contains the model parameters
#### ; Dimension of PARMOD is: DIMENSION PARMOD(10)
#### ; currently just PARMOD[0] is used for the dipole tilt angle
#### ; x,y,z input position
#### ; OUTPUT: mag .filed Bx, By, Bz at x, y, z
#### ;--------------------------------------------------------------------------------------
#### pro dipole_shielded, PARMOD, x, y, z, Bx, By, Bz
[docs]def dipole_shielded(PARMOD, x, y, z):
''' Obtain dipole field?
::
Original comments in IDL.
; WRITTEN BY K. K. KHURANA 11/2002
; MODIFIED BY H. K. SCHWARZL 11/2003
; PARMOD is an input array (real*8) that contains the model parameters
; Dimension of PARMOD is: DIMENSION PARMOD(10)
; currently just PARMOD[0] is used for the dipole tilt angle
; x,y,z input position
; OUTPUT: mag .filed Bx, By, Bz at x, y, z
'''
#### psir = PARMOD[0]
psir = PARMOD[0]
#### ;B0 = sqrt(420543d*420543d + 65920d*65920d + 24992d*24992d) ; Vip4 dipole magnitude
B0 = 426411.1411689427332021
#### B0x = B0*sin(psir)
B0x = B0*np.sin(psir)
#### B0y = 0d
B0y = 0.
#### B0z = B0*cos(psir)
B0z = B0*np.cos(psir)
####
#### ; We will first calculate the field of the dipole
#### dipole, B0x, B0y, B0z, x, y, z, Bxd, Byd, Bzd
Bxd, Byd, Bzd = dipole(B0x, B0y, B0z, x, y, z)
####
#### ; Now calculate the parallel dipole shielding field
#### if (z eq 0d) && (y eq 0d) then begin
#### cos_phi = 1d
#### sin_phi = 0d
if z == 0. and y == 0.:
cos_phi = 1.
sin_phi = 0.
else:
#### endif else begin
#### phi = atan(z,y)
#### cos_phi = cos(phi)
#### sin_phi = sin(phi)
phi = np.arctan2(z,y)
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
#### endelse
#### rho = y*cos_phi + z*sin_phi
rho = y*cos_phi + z*sin_phi
#### ; Number of dipole modes
#### Nmodes = 3
Nmodes = 3
####
#### ; Call B_mp_par(rho,phi,x,Brho1,Bphi1,Bx1,Nmodes) !no par dipole anymore
#### Brho1 = 0d
#### Bphi1 = 0d
#### Bx1 = 0d
Brho1 = 0.
Bphi1 = 0.
Bx1 = 0.
#### B_mp_perp, rho, phi, x, Brho2, Bphi2, Bx2, Nmodes
Brho2, Bphi2, Bx2 = B_mp_perp(rho, phi, x, Nmodes)
####
#### By2 = Brho2*cos_phi - Bphi2*sin_phi
#### Bz2 = Brho2*sin_phi + Bphi2*cos_phi
####
By2 = Brho2*cos_phi - Bphi2*sin_phi
Bz2 = Brho2*sin_phi + Bphi2*cos_phi
#### Bx = Bxd + Bx2
#### By = Byd + By2
#### Bz = Bzd + Bz2
Bx = Bxd + Bx2
By = Byd + By2
Bz = Bzd + Bz2
return Bx, By, Bz
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure dipole_shield_cyl_s3
#### ;--------------------------------------------------------------------------------------
#### pro dipole_shield_cyl_S3, PARMOD, ctime, rmap, pmap, zmap, Brds, Bpds, Bzds, sphiOut
[docs]def dipole_shield_cyl_S3(PARMOD, ctime, rmap, pmap, zmap, sphiOut):
#### cyl2car_pos, rmap, pmap, zmap, Xcar, Ycar, Zcar
Xcar, Ycar, Zcar = cyl2car_pos(rmap, pmap, zmap)
####
#### S32PJSO, xpJSO, ypJSO, zpJSO, Xcar, Ycar, Zcar, sphiOut
xpJSO, ypJSO, zpJSO = S32PJSO(Xcar, Ycar, Zcar, sphiOut)
#### IDL returns
#
#### dipole_shielded, PARMOD, xpJSO, ypJSO, zpJSO, BxpJSO, BypJSO, BzpJSO
(BxpJSO, BypJSO, BzpJSO) = dipole_shielded(PARMOD, xpJSO, ypJSO, zpJSO)
#### IDL returns
# PARMOD = 0. 0. 0. .... (10 elements)
# print, xpJSO, ypJSO, zpJSO, BxpJSO, BypJSO, BzpJSO
# 29.774504 -3.9380771e-08 5.7614139 8.5278990 -1.2971105e-06 -14.171122
#### Python
# 29.774497437565074, 0.0, 5.7614152250836614, 8.5279075477744062, -1.2850552970324457e-06, -14.171125158004727
#### I feel a bit worried about ypJSO==0, but otherwise, it looks ok.
_logger.debug('PARMOD = {}'.format(PARMOD))
_logger.debug('xpJSO, ypJSO, zpJSO, BxpJSO, BypJSO, BzpJSO={}'.format(str([xpJSO, ypJSO, zpJSO, BxpJSO, BypJSO, BzpJSO])))
#### PJSO2S3, BxpJSO, BypJSO, BzpJSO, BxS3, ByS3, BzS3, sphiOut
BxS3, ByS3, BzS3 = PJSO2S3(BxpJSO, BypJSO, BzpJSO, sphiOut)
####
#### car2cyl_vect, Brds, Bpds, Bzds, BxS3, ByS3, BzS3, pmap
Brds, Bpds, Bzds = car2cyl_vect(BxS3, ByS3, BzS3, pmap)
return (Brds, Bpds, Bzds)
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; funciton U
#### ;
#### ; calculates the potential terms
#### ;--------------------------------------------------------------------------------------
#### ;function U, a, c, M, rpvecP, rpvecR, sin1, sin3, cos1, sin2, M2, M2_mM
#### ; U = 0d
#### ; for i = 0, M - 1 do begin
#### ; U += total(a[i,0:M2_mM]*(rpvecP[i,M:M2])*cos1[i]*sin1[M:M2] $
#### ; + c[i,0:M2_mM]*(rpvecR[i,M:M2])*sin2[i]*sin3[M:M2])
#### ; endfor
#### ; return, U
#### ;end ; U function
####
####
#### function U_B, x1, y1, z1, a, c, p_1, p_2, r_1, r_2
####
#### Up = [0d, 0d, 0d]
#### Um = [0d, 0d, 0d]
#### B = [0d, 0d, 0d]
####
#### xp = [x1 + 0.001d, x1, x1]
#### yp = [y1, y1 + 0.001d, y1]
#### zp = [z1, z1, z1 + 0.001d]
#### xm = [x1 - 0.001d, x1, x1]
#### ym = [y1, y1 - 0.001d, y1]
#### zm = [z1, z1, z1 - 0.001d]
####
#### term_p = (term_r = dblarr(8,8))
#### for i = 0, 7 do begin
#### term_p[i,0:7] = 1d/(p_1[i]*p_1[i]) + 1d/(p_2[0:7]*p_2[0:7])
#### term_r[i,0:7] = 1d/(r_1[i]*r_1[i]) + 1d/(r_2[0:7]*r_2[0:7])
#### endfor
#### term_p = sqrt(term_p)
#### term_r = sqrt(term_r)
####
#### for i = 0, 7 do begin
#### cos_yp = cos(yp/p_1[i])
#### cos_ym = cos(ym/p_1[i])
#### sin_yp = sin(yp/r_1[i])
#### sin_ym = sin(ym/r_1[i])
#### for k = 0, 7 do begin
#### Up = Up + a[i,k]*exp(term_p[i,k]*xp)*cos_yp*sin(zp/p_2[k]) + $
#### c[i,k]*exp(term_r[i,k]*xp)*sin_yp*sin(zp/r_2[k])
####
#### Um = Um + a[i,k]*exp(term_p[i,k]*xm)*cos_ym*sin(zm/p_2[k]) + $
#### c[i,k]*exp(term_r[i,k]*xm)*sin_ym*sin(zm/r_2[k])
#### endfor
#### endfor
####
#### B = Up - Um
#### return, B
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure B_tail_shield
#### ;--------------------------------------------------------------------------------------
#### pro B_tail_shield, x, y, z, Bxout, Byout, Bzout, M, ModeIn
[docs]def B_tail_shield(x, y, z, M, ModeIn):
#### ; declare constants a, p, c, r
#### aMode = (cMode = dblarr(6, 8, 8))
aMode = np.zeros([6, 8, 8])
cMode = np.zeros([6, 8, 8])
#### p = dblarr(8, 16)
#### r = dblarr(8, 16)
p = np.zeros([8, 16])
r = np.zeros([8, 16])
####
#### ; coeff a
#### aMode[0,0,*] = [0.13217647646022754326d-20, 0.50643011679165539362d-16, -0.62311994734606299672d-13, 0.73000393989897114366d-10, $
#### -0.10249726242696675093d-07, 0.32914667786334130816d-06, -0.24800841788149776689d-05, 0.37244082145067776146d-05]
#### aMode[0,1,*] = [0.22584041437827493403d-17, 0.20229608874904410065d-12, -0.13496084238172885161d-09, 0.19512264195191708182d-06, $
#### -0.28536292479920302156d-04, 0.92636894313957363067d-03, -0.69993880953881300044d-02, 0.10518459214365802889d-01]
#### aMode[0,2,*] = [0.36523294131124144357d-14, -0.64420435736454777497d-09, -0.57096883444695674114d-07, -0.27055379124230958254d-04, $
#### 0.44272407289722757184d-02, -0.14889477801070045259d+00, 0.11356290643917379412d+01, -0.17106683807405858033d+01]
#### aMode[0,3,*] = [-0.18416420070040647516d-11, 0.27046079272174843310d-06, 0.44203411790099522704d-04, 0.12944858989375875779d-03, $
#### -0.51975317855123419619d-02, 0.20369979936316231494d+00, -0.16254799676715984801d+01, 0.24787158145629919481d+01]
#### aMode[0,4,*] = [0.13636392185761578854d-09, -0.20156479561298437097d-04, -0.39633156929923650579d-02, -0.13365044136678470021d-01, $
#### -0.10442458756051709034d-01, 0.11380270847337552453d+00, -0.88061814759885734815d+00, 0.13114182967621741404d+01]
#### aMode[0,5,*] = [-0.22305468438407642928d-08, 0.33146982079580284974d-03, 0.69006040956576439882d-01, 0.27321338618451562751d+00, $
#### 0.15256440935574380191d+00, 0.82160642450780390078d-01, -0.38978676634608163453d+00, 0.47434195602775641731d+00]
#### aMode[0,6,*] = [0.84447121731127730015d-08, -0.12569305425012970989d-02, -0.26569849211645015785d+00, -0.11069710175423441711d+01, $
#### -0.64840056779980832502d+00, -0.97562574536155111104d-02, 0.15631498731831953818d+00, -0.16203940584912954747d+00]
#### aMode[0,7,*] = [-0.63486912141303140089d-08, 0.94534756045114320954d-03, 0.20061320014222592256d+00, 0.84751706801178663397d+00, $
#### 0.50107973765980506897d+00, -0.14177955455321484379d-03, 0.56476656376811078530d-01, -0.52516865628955500255d-02]
####
#### aMode[1,0,*] = [0.13619937133704040910d-23, 0.43621985467324879692d-18, -0.15098830327683961272d-14, 0.63123957003196471404d-11, $
#### -0.11978885144275235319d-08, 0.41934992058305047279d-07, -0.31951294540338555094d-06, 0.48135741458667444803d-06]
#### aMode[1,1,*] = [0.24119037788698203250d-20, 0.40881122046673903369d-14, -0.21545235983615116381d-11, 0.11640656169328844615d-07, $
#### -0.21959806946493167778d-05, 0.76662481447507042631d-04, -0.58366072779967188566d-03, 0.87913032319456423380d-03]
#### aMode[1,2,*] = [0.31144339871300497080d-16, -0.22620785858629064435d-10, -0.74029445112059493183d-08, -0.93859788335645291112d-07, $
#### 0.82317630609657417295d-04, -0.36569184498810431982d-02, 0.29432171359264206245d-01, -0.44946922231727173269d-01]
#### aMode[1,3,*] = [-0.58080626963236774429d-14, 0.39682408099390080735d-08, 0.16840660739515575627d-05, 0.38335076286793983157d-05, $
#### 0.20836381731764990199d-03, -0.33920492396785180133d-02, 0.15820849562165701485d-01, -0.19644945238949279797d-01]
#### aMode[1,4,*] = [0.41076898294619867968d-12, -0.27989015436962660920d-06, -0.13840368413368968703d-03, -0.70869170102476930495d-03, $
#### 0.26033168308504186505d-03, -0.24526926326940952094d-02, -0.35013567157541127805d-01, 0.69275925576020052076d-01]
#### aMode[1,5,*] = [-0.70799772156343383500d-11, 0.48274807820722331896d-05, 0.24974637159216332982d-02, 0.15464887864241698700d-01, $
#### 0.11439264646130791192d-01, 0.38097770315818785036d-01, -0.13590064465091589163d+00, 0.17180997533186371128d+00]
#### aMode[1,6,*] = [0.26783558161963370025d-10, -0.18266964017400193043d-04, -0.95586630597219102156d-02, -0.62918902435418315732d-01, $
#### -0.64039968185764291064d-01, -0.22222194307812768165d-03, 0.50813426126672718297d-01, -0.46991112128761134414d-01]
#### aMode[1,7,*] = [-0.20108572992030198101d-10, 0.13715428126943252085d-04, 0.71979427421967230316d-02, 0.48199320155821636646d-01, $
#### 0.49407199622570425745d-01, 0.92225425694026483824d-02, 0.53857653772015723347d-02, 0.91367943357959795491d-03]
####
#### aMode[2,0,*] = [-0.13119415239474609968d-23, -0.60256122420913360571d-18, 0.37920759812307474057d-14, -0.18206762505356291370d-10, $
#### 0.26774304351790334521d-08, -0.87684408227198336049d-07, 0.66457031757235851543d-06, -0.99948895864718618753d-06]
#### aMode[2,1,*] = [-0.25960636282069127211d-19, -0.42000012814022671392d-16, 0.74796457397241065123d-11, -0.25318674869395683124d-07, $
#### 0.34976874196721219334d-05, -0.11242351418125537954d-03, 0.84791358984972013956d-03, -0.12736554117114038398d-02]
#### aMode[2,2,*] = [0.76343432773888704190d-16, 0.64378927275727946266d-12, 0.16598929643645348619d-08, -0.92837331289370776943d-05, $
#### 0.14285390095735186477d-02, -0.47503134695174127344d-01, 0.36153292729534660665d+00, -0.54431063875700962384d+00]
#### aMode[2,3,*] = [-0.25391252497671850107d-13, -0.31368629291300997863d-09, 0.14987071996970728449d-05, 0.17824032262389087222d-04, $
#### 0.11355427076600370650d-02, -0.37610982091768150326d-01, 0.28513697308043282063d+00, -0.42869687938445020236d+00]
#### aMode[2,4,*] = [0.17211923850683952252d-11, 0.31101556003157209140d-07, -0.12650545583107459801d-03, -0.30769692479640595728d-02, $
#### -0.45994232482610790668d-02, 0.15393025113972294448d-01, 0.86726269138905074385d-02, -0.41633915235090732664d-01]
#### aMode[2,5,*] = [-0.27638896781422706006d-10, -0.55187436783116812222d-06, 0.21199958794900246594d-02, 0.61391435739846613728d-01, $
#### 0.52329651593674135767d-01, -0.83734155027368615265d-01, 0.12671335310678004670d+00, -0.18213430508768142956d+00]
#### aMode[2,6,*] = [0.10413809392985686752d-09, 0.21324280383467479893d-05, -0.80597222736636293660d-02, -0.24311900160452664110d+00, $
#### -0.13586934681330335994d+00, 0.10699875405203342904d-02, -0.67114222283495701404d-01, 0.56104987912055399590d-02]
#### aMode[2,7,*] = [-0.78195074731598692707d-10, -0.16113432047677196834d-05, 0.60646416575606716392d-02, 0.18453527959746301334d+00, $
#### 0.97748021711249979404d-01, -0.86931770826402736673d-01, -0.25859340706079221305d-01, 0.30535054676239381521d-01]
####
#### aMode[3,0,*] = [-0.11168229337158983582d-20, -0.96888835482058599524d-16, 0.28484097342054432999d-12, -0.88991894939124396302d-09, $
#### 0.13436311859750933450d-06, -0.43949492953948432472d-05, 0.33267991354851349505d-04, -0.50016760854385964307d-04]
#### aMode[3,1,*] = [-0.20765635100251089717d-17, -0.36197386357992140659d-12, 0.35732517136988684036d-09, -0.11320141902699047964d-05, $
#### 0.16737866108943116216d-03, -0.54399134640077706492d-02, 0.41108560676656589194d-01, -0.61778272173930721677d-01]
#### aMode[3,2,*] = [-0.49227064641270947831d-14, 0.13773911494082902162d-08, 0.28000536764000152345d-06, -0.60982707102378110874d-05, $
#### -0.45545917365779216012d-03, 0.30474617144202409413d-01, -0.26240435518843665541d+00, 0.40666851745242809101d+00]
#### aMode[3,3,*] = [0.24046710130055206633d-11, -0.56928231643345874601d-06, -0.14658425595003152785d-03, -0.54669311790752939117d-03, $
#### 0.63289641147589756897d-01, -0.23416266019622873351d+01, 0.18376229679993183907d+02, -0.27890706711124808592d+02]
#### aMode[3,4,*] = [-0.17790734879856566763d-09, 0.41757078162568088686d-04, 0.12661549470502249103d-01, 0.30791263485754920559d-01, $
#### 0.52395704373591174274d-01, -0.15856752943903087427d+01, 0.11848914793373501730d+02, -0.17521753728113509396d+02]
#### aMode[3,5,*] = [0.29086256117761690731d-08, -0.68275932316891383422d-03, -0.21784801486997875663d+00, -0.66623989758662736093d+00, $
#### -0.28050159504711622560d+00, -0.67774272696429793683d+00, 0.10271422822879689995d+01, -0.63392178592402839143d+00]
#### aMode[3,6,*] = [-0.11009990780568283952d-07, 0.25849182195501789749d-02, 0.83614032826616710991d+00, 0.27685938326621157834d+01, $
#### 0.16426046025686344975d+01, 0.78572426172170208857d+00, -0.18258375964853595263d+01, 0.12152196449017165225d+01]
#### aMode[3,7,*] = [0.82768727273709110647d-08, -0.19433482187601436308d-02, -0.63081127374393464180d+00, -0.21366815638955616307d+01, $
#### -0.13446774521063316054d+01, -0.58592309940003861612d+00, 0.31827661237739279798d+00, -0.20240366730373042791d+00]
####
#### aMode[4,0,*] = [-0.13947208148057608312d-20, -0.20210120268149411870d-15, 0.45539475964762736737d-12, -0.11170239869341296312d-08, $
#### 0.16610385153669340319d-06, -0.54155353001624000341d-05, 0.40961114502502580236d-04, -0.61571027127369957199d-04]
#### aMode[4,1,*] = [-0.15010983747727246750d-17, -0.56103448860629763217d-12, 0.43831990483680147718d-09, -0.10609291877688533656d-05, $
#### 0.15387190085182278487d-03, -0.49799395553291745386d-02, 0.37594128280158018995d-01, -0.56482482158755402679d-01]
#### aMode[4,2,*] = [-0.84949700310157272298d-14, 0.19131227654780431635d-08, 0.23797873673175664599d-06, 0.67357172013098089990d-04, $
#### -0.11412327840147937774d-01, 0.38787752218711162299d+00, -0.29661484451688613361d+01, 0.44710237927891309794d+01]
#### aMode[4,3,*] = [0.36903300990875096410d-11, -0.78266544732189755606d-06, -0.14356129913638142170d-03, -0.52616304038202654780d-03, $
#### 0.50792590874790795041d-01, -0.18789041423837609556d+01, 0.14745664699188227864d+02, -0.22381820221529649117d+02]
#### aMode[4,4,*] = [-0.26744771000491356360d-09, 0.57361959249991540943d-04, 0.12466277937525278574d-01, 0.35752121631968263315d-01, $
#### 0.49443768069258187125d-01, -0.11592738471082895124d+01, 0.85575941704881230975d+01, -0.12621255590568500881d+02]
#### aMode[4,5,*] = [0.43525625631567042006d-08, -0.93805114409436534117d-03, -0.21484367014734804257d+00, -0.75606598389073855770d+00, $
#### -0.38916312244605659742d+00, -0.45744351551847461934d+00, 0.91714655267490634571d+00, -0.70503052173487583687d+00]
#### aMode[4,6,*] = [-0.16457579884296458239d-07, 0.35517012004150618764d-02, 0.82496531512910689087d+00, 0.31057726963483678339d+01, $
#### 0.18813656203693941648d+01, 0.51092359576032500001d+00, -0.12428990388842335867d+01, 0.89617638629098603786d+00]
#### aMode[4,7,*] = [0.12368783158269132105d-07, -0.26702315450896003667d-02, -0.62244835948407777337d+00, -0.23879836056356684714d+01, $
#### -0.14962346641346604414d+01, -0.32180960449888429408d+00, 0.14080714553385571541d+00, -0.12793599129792307955d+00]
####
#### aMode[5,0,*] = [-0.71007031654422592126d-20, -0.82331351329343824829d-15, 0.14860700245019204501d-11, -0.18336776863420102046d-08, $
#### 0.26381273327308489839d-06, -0.85305767261458509409d-05, 0.64390806859274407614d-04, -0.96740352003902945199d-04]
#### aMode[5,1,*] = [-0.47759968145964668551d-17, -0.11110703618965760419d-11, 0.93331659872542349631d-09, -0.11361217227825839426d-05, $
#### 0.16101720194965896126d-03, -0.51853068878838977084d-02, 0.39098877702660601585d-01, -0.58726418277052969685d-01]
#### aMode[5,2,*] = [-0.20561632805425666958d-13, 0.30615328304455808883d-08, 0.24271749535760407390d-06, 0.15937085189299786236d-03, $
#### -0.25318551119527441528d-01, 0.84552595243637060917d+00, -0.64376874089793245659d+01, 0.96932895116823303283d+01]
#### aMode[5,3,*] = [0.87503551825887093684d-11, -0.13163374744119336057d-05, -0.21315095274792228430d-03, -0.75679626156015320503d-03, $
#### 0.43875888066817978483d-01, -0.16187067723465442981d+01, 0.12704254977048992092d+02, -0.19288846975099222191d+02]
#### aMode[5,4,*] = [-0.62980025164513007140d-09, 0.99033147185834717873d-04, 0.19286483194829802556d-01, 0.67190391627066974322d-01, $
#### 0.59018569821895905391d-01, -0.57675212643809139478d+00, 0.41759144712554761014d+01, -0.60616007357605896643d+01]
#### aMode[5,5,*] = [0.10235134109675949609d-07, -0.16329320547735449054d-02, -0.33676325501000334838d+00, -0.13741116541224409619d+01, $
#### -0.82897929864310349046d+00, -0.29064544069947594095d+00, 0.99716723490406433683d+00, -0.11955304085306053352d+01]
#### aMode[5,6,*] = [-0.38688354152446052580d-07, 0.61964020776766650655d-02, 0.12977027325262353585d+01, 0.55698639583852660450d+01, $
#### 0.35055306218486181890d+01, 0.24501453852569592406d+00, -0.51229945875859312920d+00, 0.42642747170716877036d+00]
#### aMode[5,7,*] = [0.29074290348871754119d-07, -0.46611910769336493132d-02, -0.98002188003461245813d+00, -0.42652335402388423801d+01, $
#### -0.27213526220765031915d+01, -0.21604814310942872523d+00, -0.40325674485435270000d+00, 0.12987860657487151350d+00]
####
aMode[0,0,:] = [0.13217647646022754326e-20, 0.50643011679165539362e-16, -0.62311994734606299672e-13, 0.73000393989897114366e-10,
-0.10249726242696675093e-07, 0.32914667786334130816e-06, -0.24800841788149776689e-05, 0.37244082145067776146e-05]
aMode[0,1,:] = [0.22584041437827493403e-17, 0.20229608874904410065e-12, -0.13496084238172885161e-09, 0.19512264195191708182e-06,
-0.28536292479920302156e-04, 0.92636894313957363067e-03, -0.69993880953881300044e-02, 0.10518459214365802889e-01]
aMode[0,2,:] = [0.36523294131124144357e-14, -0.64420435736454777497e-09, -0.57096883444695674114e-07, -0.27055379124230958254e-04,
0.44272407289722757184e-02, -0.14889477801070045259e+00, 0.11356290643917379412e+01, -0.17106683807405858033e+01]
aMode[0,3,:] = [-0.18416420070040647516e-11, 0.27046079272174843310e-06, 0.44203411790099522704e-04, 0.12944858989375875779e-03,
-0.51975317855123419619e-02, 0.20369979936316231494e+00, -0.16254799676715984801e+01, 0.24787158145629919481e+01]
aMode[0,4,:] = [0.13636392185761578854e-09, -0.20156479561298437097e-04, -0.39633156929923650579e-02, -0.13365044136678470021e-01,
-0.10442458756051709034e-01, 0.11380270847337552453e+00, -0.88061814759885734815e+00, 0.13114182967621741404e+01]
aMode[0,5,:] = [-0.22305468438407642928e-08, 0.33146982079580284974e-03, 0.69006040956576439882e-01, 0.27321338618451562751e+00,
0.15256440935574380191e+00, 0.82160642450780390078e-01, -0.38978676634608163453e+00, 0.47434195602775641731e+00]
aMode[0,6,:] = [0.84447121731127730015e-08, -0.12569305425012970989e-02, -0.26569849211645015785e+00, -0.11069710175423441711e+01,
-0.64840056779980832502e+00, -0.97562574536155111104e-02, 0.15631498731831953818e+00, -0.16203940584912954747e+00]
aMode[0,7,:] = [-0.63486912141303140089e-08, 0.94534756045114320954e-03, 0.20061320014222592256e+00, 0.84751706801178663397e+00,
0.50107973765980506897e+00, -0.14177955455321484379e-03, 0.56476656376811078530e-01, -0.52516865628955500255e-02]
aMode[1,0,:] = [0.13619937133704040910e-23, 0.43621985467324879692e-18, -0.15098830327683961272e-14, 0.63123957003196471404e-11,
-0.11978885144275235319e-08, 0.41934992058305047279e-07, -0.31951294540338555094e-06, 0.48135741458667444803e-06]
aMode[1,1,:] = [0.24119037788698203250e-20, 0.40881122046673903369e-14, -0.21545235983615116381e-11, 0.11640656169328844615e-07,
-0.21959806946493167778e-05, 0.76662481447507042631e-04, -0.58366072779967188566e-03, 0.87913032319456423380e-03]
aMode[1,2,:] = [0.31144339871300497080e-16, -0.22620785858629064435e-10, -0.74029445112059493183e-08, -0.93859788335645291112e-07,
0.82317630609657417295e-04, -0.36569184498810431982e-02, 0.29432171359264206245e-01, -0.44946922231727173269e-01]
aMode[1,3,:] = [-0.58080626963236774429e-14, 0.39682408099390080735e-08, 0.16840660739515575627e-05, 0.38335076286793983157e-05,
0.20836381731764990199e-03, -0.33920492396785180133e-02, 0.15820849562165701485e-01, -0.19644945238949279797e-01]
aMode[1,4,:] = [0.41076898294619867968e-12, -0.27989015436962660920e-06, -0.13840368413368968703e-03, -0.70869170102476930495e-03,
0.26033168308504186505e-03, -0.24526926326940952094e-02, -0.35013567157541127805e-01, 0.69275925576020052076e-01]
aMode[1,5,:] = [-0.70799772156343383500e-11, 0.48274807820722331896e-05, 0.24974637159216332982e-02, 0.15464887864241698700e-01,
0.11439264646130791192e-01, 0.38097770315818785036e-01, -0.13590064465091589163e+00, 0.17180997533186371128e+00]
aMode[1,6,:] = [0.26783558161963370025e-10, -0.18266964017400193043e-04, -0.95586630597219102156e-02, -0.62918902435418315732e-01,
-0.64039968185764291064e-01, -0.22222194307812768165e-03, 0.50813426126672718297e-01, -0.46991112128761134414e-01]
aMode[1,7,:] = [-0.20108572992030198101e-10, 0.13715428126943252085e-04, 0.71979427421967230316e-02, 0.48199320155821636646e-01,
0.49407199622570425745e-01, 0.92225425694026483824e-02, 0.53857653772015723347e-02, 0.91367943357959795491e-03]
aMode[2,0,:] = [-0.13119415239474609968e-23, -0.60256122420913360571e-18, 0.37920759812307474057e-14, -0.18206762505356291370e-10,
0.26774304351790334521e-08, -0.87684408227198336049e-07, 0.66457031757235851543e-06, -0.99948895864718618753e-06]
aMode[2,1,:] = [-0.25960636282069127211e-19, -0.42000012814022671392e-16, 0.74796457397241065123e-11, -0.25318674869395683124e-07,
0.34976874196721219334e-05, -0.11242351418125537954e-03, 0.84791358984972013956e-03, -0.12736554117114038398e-02]
aMode[2,2,:] = [0.76343432773888704190e-16, 0.64378927275727946266e-12, 0.16598929643645348619e-08, -0.92837331289370776943e-05,
0.14285390095735186477e-02, -0.47503134695174127344e-01, 0.36153292729534660665e+00, -0.54431063875700962384e+00]
aMode[2,3,:] = [-0.25391252497671850107e-13, -0.31368629291300997863e-09, 0.14987071996970728449e-05, 0.17824032262389087222e-04,
0.11355427076600370650e-02, -0.37610982091768150326e-01, 0.28513697308043282063e+00, -0.42869687938445020236e+00]
aMode[2,4,:] = [0.17211923850683952252e-11, 0.31101556003157209140e-07, -0.12650545583107459801e-03, -0.30769692479640595728e-02,
-0.45994232482610790668e-02, 0.15393025113972294448e-01, 0.86726269138905074385e-02, -0.41633915235090732664e-01]
aMode[2,5,:] = [-0.27638896781422706006e-10, -0.55187436783116812222e-06, 0.21199958794900246594e-02, 0.61391435739846613728e-01,
0.52329651593674135767e-01, -0.83734155027368615265e-01, 0.12671335310678004670e+00, -0.18213430508768142956e+00]
aMode[2,6,:] = [0.10413809392985686752e-09, 0.21324280383467479893e-05, -0.80597222736636293660e-02, -0.24311900160452664110e+00,
-0.13586934681330335994e+00, 0.10699875405203342904e-02, -0.67114222283495701404e-01, 0.56104987912055399590e-02]
aMode[2,7,:] = [-0.78195074731598692707e-10, -0.16113432047677196834e-05, 0.60646416575606716392e-02, 0.18453527959746301334e+00,
0.97748021711249979404e-01, -0.86931770826402736673e-01, -0.25859340706079221305e-01, 0.30535054676239381521e-01]
aMode[3,0,:] = [-0.11168229337158983582e-20, -0.96888835482058599524e-16, 0.28484097342054432999e-12, -0.88991894939124396302e-09,
0.13436311859750933450e-06, -0.43949492953948432472e-05, 0.33267991354851349505e-04, -0.50016760854385964307e-04]
aMode[3,1,:] = [-0.20765635100251089717e-17, -0.36197386357992140659e-12, 0.35732517136988684036e-09, -0.11320141902699047964e-05,
0.16737866108943116216e-03, -0.54399134640077706492e-02, 0.41108560676656589194e-01, -0.61778272173930721677e-01]
aMode[3,2,:] = [-0.49227064641270947831e-14, 0.13773911494082902162e-08, 0.28000536764000152345e-06, -0.60982707102378110874e-05,
-0.45545917365779216012e-03, 0.30474617144202409413e-01, -0.26240435518843665541e+00, 0.40666851745242809101e+00]
aMode[3,3,:] = [0.24046710130055206633e-11, -0.56928231643345874601e-06, -0.14658425595003152785e-03, -0.54669311790752939117e-03,
0.63289641147589756897e-01, -0.23416266019622873351e+01, 0.18376229679993183907e+02, -0.27890706711124808592e+02]
aMode[3,4,:] = [-0.17790734879856566763e-09, 0.41757078162568088686e-04, 0.12661549470502249103e-01, 0.30791263485754920559e-01,
0.52395704373591174274e-01, -0.15856752943903087427e+01, 0.11848914793373501730e+02, -0.17521753728113509396e+02]
aMode[3,5,:] = [0.29086256117761690731e-08, -0.68275932316891383422e-03, -0.21784801486997875663e+00, -0.66623989758662736093e+00,
-0.28050159504711622560e+00, -0.67774272696429793683e+00, 0.10271422822879689995e+01, -0.63392178592402839143e+00]
aMode[3,6,:] = [-0.11009990780568283952e-07, 0.25849182195501789749e-02, 0.83614032826616710991e+00, 0.27685938326621157834e+01,
0.16426046025686344975e+01, 0.78572426172170208857e+00, -0.18258375964853595263e+01, 0.12152196449017165225e+01]
aMode[3,7,:] = [0.82768727273709110647e-08, -0.19433482187601436308e-02, -0.63081127374393464180e+00, -0.21366815638955616307e+01,
-0.13446774521063316054e+01, -0.58592309940003861612e+00, 0.31827661237739279798e+00, -0.20240366730373042791e+00]
aMode[4,0,:] = [-0.13947208148057608312e-20, -0.20210120268149411870e-15, 0.45539475964762736737e-12, -0.11170239869341296312e-08,
0.16610385153669340319e-06, -0.54155353001624000341e-05, 0.40961114502502580236e-04, -0.61571027127369957199e-04]
aMode[4,1,:] = [-0.15010983747727246750e-17, -0.56103448860629763217e-12, 0.43831990483680147718e-09, -0.10609291877688533656e-05,
0.15387190085182278487e-03, -0.49799395553291745386e-02, 0.37594128280158018995e-01, -0.56482482158755402679e-01]
aMode[4,2,:] = [-0.84949700310157272298e-14, 0.19131227654780431635e-08, 0.23797873673175664599e-06, 0.67357172013098089990e-04,
-0.11412327840147937774e-01, 0.38787752218711162299e+00, -0.29661484451688613361e+01, 0.44710237927891309794e+01]
aMode[4,3,:] = [0.36903300990875096410e-11, -0.78266544732189755606e-06, -0.14356129913638142170e-03, -0.52616304038202654780e-03,
0.50792590874790795041e-01, -0.18789041423837609556e+01, 0.14745664699188227864e+02, -0.22381820221529649117e+02]
aMode[4,4,:] = [-0.26744771000491356360e-09, 0.57361959249991540943e-04, 0.12466277937525278574e-01, 0.35752121631968263315e-01,
0.49443768069258187125e-01, -0.11592738471082895124e+01, 0.85575941704881230975e+01, -0.12621255590568500881e+02]
aMode[4,5,:] = [0.43525625631567042006e-08, -0.93805114409436534117e-03, -0.21484367014734804257e+00, -0.75606598389073855770e+00,
-0.38916312244605659742e+00, -0.45744351551847461934e+00, 0.91714655267490634571e+00, -0.70503052173487583687e+00]
aMode[4,6,:] = [-0.16457579884296458239e-07, 0.35517012004150618764e-02, 0.82496531512910689087e+00, 0.31057726963483678339e+01,
0.18813656203693941648e+01, 0.51092359576032500001e+00, -0.12428990388842335867e+01, 0.89617638629098603786e+00]
aMode[4,7,:] = [0.12368783158269132105e-07, -0.26702315450896003667e-02, -0.62244835948407777337e+00, -0.23879836056356684714e+01,
-0.14962346641346604414e+01, -0.32180960449888429408e+00, 0.14080714553385571541e+00, -0.12793599129792307955e+00]
aMode[5,0,:] = [-0.71007031654422592126e-20, -0.82331351329343824829e-15, 0.14860700245019204501e-11, -0.18336776863420102046e-08,
0.26381273327308489839e-06, -0.85305767261458509409e-05, 0.64390806859274407614e-04, -0.96740352003902945199e-04]
aMode[5,1,:] = [-0.47759968145964668551e-17, -0.11110703618965760419e-11, 0.93331659872542349631e-09, -0.11361217227825839426e-05,
0.16101720194965896126e-03, -0.51853068878838977084e-02, 0.39098877702660601585e-01, -0.58726418277052969685e-01]
aMode[5,2,:] = [-0.20561632805425666958e-13, 0.30615328304455808883e-08, 0.24271749535760407390e-06, 0.15937085189299786236e-03,
-0.25318551119527441528e-01, 0.84552595243637060917e+00, -0.64376874089793245659e+01, 0.96932895116823303283e+01]
aMode[5,3,:] = [0.87503551825887093684e-11, -0.13163374744119336057e-05, -0.21315095274792228430e-03, -0.75679626156015320503e-03,
0.43875888066817978483e-01, -0.16187067723465442981e+01, 0.12704254977048992092e+02, -0.19288846975099222191e+02]
aMode[5,4,:] = [-0.62980025164513007140e-09, 0.99033147185834717873e-04, 0.19286483194829802556e-01, 0.67190391627066974322e-01,
0.59018569821895905391e-01, -0.57675212643809139478e+00, 0.41759144712554761014e+01, -0.60616007357605896643e+01]
aMode[5,5,:] = [0.10235134109675949609e-07, -0.16329320547735449054e-02, -0.33676325501000334838e+00, -0.13741116541224409619e+01,
-0.82897929864310349046e+00, -0.29064544069947594095e+00, 0.99716723490406433683e+00, -0.11955304085306053352e+01]
aMode[5,6,:] = [-0.38688354152446052580e-07, 0.61964020776766650655e-02, 0.12977027325262353585e+01, 0.55698639583852660450e+01,
0.35055306218486181890e+01, 0.24501453852569592406e+00, -0.51229945875859312920e+00, 0.42642747170716877036e+00]
aMode[5,7,:] = [0.29074290348871754119e-07, -0.46611910769336493132e-02, -0.98002188003461245813e+00, -0.42652335402388423801e+01,
-0.27213526220765031915e+01, -0.21604814310942872523e+00, -0.40325674485435270000e+00, 0.12987860657487151350e+00]
####
#### ; coeff c
#### cMode[0,0,*] = [0.62608370709224905326d-18, -0.69296927042311144973d-15, 0.53793142672880289723d-22, 0.55953376048336647130d-13, $
#### -0.13644845214970875435d-11, 0.51560868494277940499d-10, -0.38841223167697958018d-09, 0.58169693048278663383d-09]
#### cMode[0,1,*] = [0.79948744245788159190d-02, 0.28880230453715999061d-01, 0.18816017618889148366d-03, 0.30432220614412108794d-01, $
#### 0.20335251312180755434d-01, 0.12253559136086547010d+01, -0.11162104253917091156d+02, 0.17710281832472645646d+02]
#### cMode[0,2,*] = [-0.30562966462902028119d+03, -0.11990795349370895195d+04, -0.65802994518949660118d+01, -0.12273805927634358070d+04, $
#### -0.65349097657709229736d+03, -0.12824769943213845024d+03, -0.88590404530202224719d+02, -0.92014824378423334394d+02]
#### cMode[0,3,*] = [-0.95152195796331646704d+02, -0.36956353969242154988d+03, -0.20695307247843501841d+01, -0.38362497103275590148d+03, $
#### -0.22776414558915512031d+03, -0.46624811555373160132d+02, 0.45685596890036768158d+02, -0.25467538299896208542d+03]
#### cMode[0,4,*] = [0.32918009134216310584d+03, 0.12851583684685703445d+04, 0.71222223072393315845d+01, 0.13249069943127731452d+04, $
#### 0.74436590714498462872d+03, 0.14749387222546454623d+03, 0.42968129090712094964d+02, 0.33666974003731491293d+03]
#### cMode[0,5,*] = [0.20156315877857942098d+03, 0.80514059433962259504d+03, 0.42634222620010087112d+01, 0.79880139299917756190d+03, $
#### 0.34272486825284445011d+03, 0.52606372363987690121d+02, 0.74826789775518696146d+01, -0.49206154796489274261d+01]
#### cMode[0,6,*] = [-0.53421252140832686805d+03, -0.21573846013695305856d+04, -0.11180286753450849879d+02, -0.20919394990520148169d+04, $
#### -0.77083430972020012816d+03, -0.74269446238765679524d+02, 0.86081879160377337001d+01, -0.33189293006747679903d+01]
#### cMode[0,7,*] = [0.64042667076955872573d+03, 0.25941583461638844099d+04, 0.13364255832056319839d+02, 0.24982768986941823463d+04, $
#### 0.87945971890665362025d+03, 0.70070351276938351858d+02, -0.86139878258589916981d+01, 0.30011846979293261838d+01]
####
#### cMode[1,0,*] = [-0.63013334128716040893d-13, 0.60865772431920603935d-05, -0.40016262502772708131d-09, -0.23069864565108759713d-04, $
#### 0.23255146051530801720d-03,-0.51951144258714805346d-02, 0.36136455088219365805d-01, -0.53017682607232261560d-01]
#### cMode[1,1,*] = [0.55436122531818750047d-03, -0.66992043381081387565d-02, 0.83803970952918511727d-02, 0.80754210845424374554d-02, $
#### -0.61427751801047527635d-02, 0.54122156026262873140d+00, -0.43010880044939874267d+01, 0.65777678112627633311d+01]
#### cMode[1,2,*] = [-0.23084376754368056694d+02, 0.20899246010739616075d+03, -0.36836390749521057408d+03, -0.52914996781002994197d+03, $
#### -0.20007855689466294002d+03, -0.36056061195965001253d+02, -0.16132602175333610183d+02, 0.21643026623710106548d+01]
#### cMode[1,3,*] = [-0.30489331973819799870d+01, 0.33304902492752548326d+02, -0.47908018506183607243d+02, -0.66968634395671120529d+02, $
#### -0.34606055921065168590d+02, -0.74942275618948936966d+01, 0.20818848022147209420d+02, -0.68532587218682348151d+02]
#### cMode[1,4,*] = [0.44092752992499013586d+01, -0.47413797913825259655d+02, 0.69401964003076273002d+02, 0.97614272148050673649d+02, $
#### 0.48704912842393142113d+02, 0.99941378554753921292d+01, -0.14570861805562320689d+02, 0.65068431858518227528d+02]
#### cMode[1,5,*] = [0.27754841758044097588d+02, -0.24286715291266780525d+03, 0.44379959765035561503d+03, 0.63697225078825798760d+03, $
#### 0.23079424266757904149d+03, 0.39146925311929221535d+02, 0.12991883575647105075d+02, -0.34552170740359109402d+01]
#### cMode[1,6,*] = [-0.20272909693943255149d+02, 0.14757207634388509021d+03, -0.32708628601917539846d+03, -0.46225115463536488036d+03, $
#### -0.13753186496272293837d+03, -0.14358017591372269627d+02, 0.41733787941502891172d+00, -0.14829734668187728452d+00]
#### cMode[1,7,*] = [0.22054195737479753702d+02, -0.15104824369566767217d+03, 0.35667767962821925742d+03, 0.50036616092734789162d+03, $
#### 0.14078440554104824755d+03, 0.12025910667067751802d+02, -0.70328949883151059552d+00, 0.24434145951697048282d+00]
####
#### cMode[2,0,*] = [-0.28862183650305684778d-01, 0.16202748718849406373d-02, 0.62452442623701900359d-08, -0.10292552636752272388d-04, $
#### 0.14899979598450072693d-10, 0.54701873831866327790d+00, -0.33801096158365382393d+01, 0.46663975792709919687d+01]
#### cMode[2,1,*] = [0.82086900572783427776d+00, 0.66781995371601139410d+00, 0.17689517609606757453d+01, 0.18458732991639696052d+01, $
#### -0.13408206456871818446d-02, -0.55045820392109447993d+01, 0.34057914707574474810d+02, -0.44172744618314503384d+02]
#### cMode[2,2,*] = [-0.33440561402073085695d+04, -0.63006668870531044035d+04, -0.16148889754716360123d+05, -0.16734685574071548330d+05, $
#### 0.98894432761575998824d+01, -0.10391741179165414621d+04, -0.17857836173455567951d+04, 0.29803489709704852117d+04]
#### cMode[2,3,*] = [0.31915064954348162373d+04, 0.59270193788343030760d+04, 0.15182397380443129364d+05, 0.15734311019448650625d+05, $
#### -0.93434105903445914265d+01, 0.99260524095183413351d+03, 0.20341329769010583206d+04, -0.32962143597458424260d+04]
#### cMode[2,4,*] = [-0.92169284592923155230d+02, -0.14814666530740792538d+03, -0.37839662711403616590d+03, -0.39255805350361141492d+03, $
#### 0.24567231705288499199d+00, -0.11195431933835682247d+02, -0.27154689478860825069d+03, 0.37730991124653248114d+03]
#### cMode[2,5,*] = [0.54075863200485683179d+03, 0.13068947201920562140d+04, 0.34015030753731947399d+04, 0.35222027482034166112d+04, $
#### -0.19372427035775936943d+01, 0.10115813738135293053d+03, -0.23832476149210184424d+02, -0.16890808415607296844d+02]
#### cMode[2,6,*] = [-0.77613909679256050111d+03, -0.22744934351169181496d+04, -0.60171788754328519033d+04, -0.62294563978631467549d+04, $
#### 0.32775011290202917813d+01, -0.72131812461142832404d+02, 0.27116603294455088324d+02, -0.71224576605140779150d+01]
#### cMode[2,7,*] = [0.69369255218888747904d+03, 0.22078310233442990373d+04, 0.58858136079721381506d+04, 0.60932655114670426499d+04, $
#### -0.31519955877359868701d+01, 0.44446156981790601037d+02, -0.14233922505288183479d+02, 0.34387488557126792976d+01]
####
#### cMode[3,0,*] = [-0.36868002035455376130d-21, 0.10541748946585167701d-17, 0.27755566547052099579d-13, -0.34780372211895671519d-15, $
#### -0.99328472275559516191d-12, 0.37378357885053432596d-10, -0.27315035342737905565d-09, 0.40965931239852784173d-09]
#### cMode[3,1,*] = [-0.33737087677076624814d-04, -0.11396625664731094840d-02, -0.10789923154822660400d-01, -0.80796849963401182748d-02, $
#### 0.32166227895723737972d-01, -0.25571424523548520468d+01, 0.21267119641439142796d+02, -0.33078925310689939465d+02]
#### cMode[3,2,*] = [-0.17479871537929662750d+01, -0.72444600080193755076d+02, -0.55033093275560522883d+03, -0.50210145669063086515d+03, $
#### -0.25844521786360195036d+03, -0.55113408951058726614d+02, 0.20187135793924890769d+03, -0.81289982453232934034d+03]
#### cMode[3,3,*] = [0.78603614447709357904d+00, 0.32329428441630612134d+02, 0.24762820448825531016d+03, 0.22413360241783610860d+03, $
#### 0.11996249854833778059d+03, 0.30448969122459486058d+02, -0.23988687682537666034d+03, 0.67671098556149500424d+03]
#### cMode[3,4,*] = [0.12159934205840530196d+01, 0.51016427583112839982d+02, 0.38232595903016939331d+03, 0.35346430190826185757d+03, $
#### 0.17043627550273194870d+03, 0.36343814970370997841d+02, 0.25302217119464871508d+02, 0.16017628158741096910d+03]
#### cMode[3,5,*] = [-0.13415948717761450037d+01, -0.59032589809690589888d+02, -0.41777999557453684431d+03, -0.40902457469604014406d+03, $
#### -0.14475183395103458749d+03, -0.29478536348031818548d+02, -0.13387094046024579085d+01, 0.23831368497822609242d+01]
#### cMode[3,6,*] = [0.55621147775238517496d+01, 0.24890779339281028370d+03, 0.17228008107309571883d+04, 0.17256605312262474072d+04, $
#### 0.52806415548478540245d+03, 0.64772764588200661961d+02, -0.12923255027738644784d+02, 0.58752559976531326668d+01]
#### cMode[3,7,*] = [-0.73285857760141963623d+01, -0.32941251542733325230d+03, -0.22662989686584209536d+04, -0.22842752627361360140d+04, $
#### -0.66988125859857134969d+03, -0.67321294156948754405d+02, 0.12094082762406057618d+02, -0.50344380663030463551d+01]
####
#### cMode[4,0,*] = [0.57951308146739251014d+01, -0.14999087019121093433d+03, -0.19559326802070165385d+01, -0.66720804039446477418d+03, $
#### -0.75756207610197607849d+02, -0.79695623597154439110d+02, 0.12741669724908954997d+04, -0.25779398744795645193d+04]
#### cMode[4,1,*] = [0.25559373539070939784d+03, -0.26954324048080380293d+04, 0.42741834334120607508d+02, -0.13609107834915077361d+05, $
#### -0.18219129315199742435d+04, -0.49540503909243787106d+03, -0.86525145704320394202d+03, -0.19349682251463748983d+04]
#### cMode[4,2,*] = [-0.87923543175638503299d+05, 0.30861003822923382955d+06, -0.25217316373834348652d+05, 0.21333098058144601694d+07, $
#### 0.33404584481156045505d+06, 0.42354436002087458845d+05, 0.58637814837743640339d+04, -0.30291614862820930298d+04]
#### cMode[4,3,*] = [0.21501887419077996277d+06, -0.98416860975998154970d+06, 0.50756266605703403982d+05, -0.58885368624678218196d+07, $
#### -0.88230632458050202160d+06, -0.24671074506661172520d+06, -0.14496629047056444505d+06, 0.68193611837692298394d+05]
#### cMode[4,4,*] = [0.18016800297232990146d+02, 0.11854731605905291047d+00, 0.41242056558479642802d-02, 0.10788775915866732901d+00, $
#### 0.76539029541206167195d+00, 0.67292188204585867694d+02, -0.31151519702859546967d+03, 0.44202304203357769551d+03]
#### cMode[4,5,*] = [0.20068140169309098830d+06, -0.93099363016509286694d+06, 0.47215756436253171202d+05, -0.55473705494319487385d+07, $
#### -0.82963829864856943885d+06, -0.23595446376207669381d+06, -0.14387731979865256981d+06, 0.66726697732453104094d+05]
#### cMode[4,6,*] = [-0.10659833147754620430d+06, 0.36525155976544767533d+06, -0.32972599682804299980d+05, 0.26280915665374022793d+07, $
#### 0.41399542203693986408d+06, 0.40200322950430580348d+05, 0.33714384362519780324d+04, -0.16170695162643912823d+04]
#### cMode[4,7,*] = [0.45385852289903034773d+05, -0.15367285853252436567d+06, 0.15938904084451054998d+05, -0.11759967065602889846d+07, $
#### -0.18617563074881875451d+06, -0.11629254312624910383d+05, -0.37356632983270783299d+03, 0.13368622059973152005d+03]
####
#### cMode[5,0,*] = [0.37071935082976326114d-15, 0.46201169190899653571d-08, -0.62075147292222574435d-11, -0.36146708876274598054d-06, $
#### 0.87046861427880060091d-05, -0.10076975633470059978d-03, 0.89298240723925399464d-03, -0.14060364310276510124d-02]
#### cMode[5,1,*] = [-0.56188781974792609830d-01, -0.40650854381275376425d+01, -0.14988738284605299000d+01, -0.38632728626301734209d+01, $
#### -0.22261011246466999580d+01, -0.14212548821848476343d+02, 0.19383430976405556123d+03, -0.34806064213022160913d+03]
#### cMode[5,2,*] = [0.27593841470781566016d+03, 0.22140385471263797079d+05, 0.78362796461845718454d+04, 0.19668830633326862766d+05, $
#### 0.10181292569179314355d+05, 0.39096643822536609746d+04, 0.26876686999693699675d+04, 0.33310618692009836827d+04]
#### cMode[5,3,*] = [0.44061551736067183782d+02, 0.34746817101944778016d+04, 0.12384774155446498511d+04, 0.31312056838986070950d+04, $
#### 0.17260274467090328087d+04, 0.70922681375271823256d+03, -0.30968222241454737009d+03, 0.29976750308301927105d+04]
#### cMode[5,4,*] = [-0.29826526924997253331d+03, -0.23828019858653441964d+05, -0.84485266842828430355d+04, -0.21247202725657969857d+05, $
#### -0.11181459854792739072d+05, -0.43708622712852429614d+04, -0.24560179552632419586d+04, -0.60645360750106398484d+04]
#### cMode[5,5,*] = [-0.74893303991611572811d+02, -0.62909773314430736945d+04, -0.21848887114479706994d+04, -0.53368828632132494504d+04, $
#### -0.22162146435679854761d+04, -0.52579256909996399116d+03, -0.11309723650631768876d+03, 0.63863048387678587047d+02]
#### cMode[5,6,*] = [0.19717911116782794067d+03, 0.16899744672548959734d+05, 0.58199397638247170050d+04, 0.13983841291452328015d+05, $
#### 0.50676093636595025415d+04, 0.75432123378980122652d+03, 0.10996364977646411276d+02, -0.48990730626876937137d+01]
#### cMode[5,7,*] = [-0.22796489432205970793d+03, -0.19665709629971770411d+05, -0.67538957248075544015d+04, -0.16130393271007925193d+05, $
#### -0.55557948391593710013d+04, -0.67829337377268315023d+03, 0.63691256079676055179d+01, -0.36820923261738309761d+01]
cMode[0,0,:] = [0.62608370709224905326e-18, -0.69296927042311144973e-15, 0.53793142672880289723e-22, 0.55953376048336647130e-13,
-0.13644845214970875435e-11, 0.51560868494277940499e-10, -0.38841223167697958018e-09, 0.58169693048278663383e-09]
cMode[0,1,:] = [0.79948744245788159190e-02, 0.28880230453715999061e-01, 0.18816017618889148366e-03, 0.30432220614412108794e-01,
0.20335251312180755434e-01, 0.12253559136086547010e+01, -0.11162104253917091156e+02, 0.17710281832472645646e+02]
cMode[0,2,:] = [-0.30562966462902028119e+03, -0.11990795349370895195e+04, -0.65802994518949660118e+01, -0.12273805927634358070e+04,
-0.65349097657709229736e+03, -0.12824769943213845024e+03, -0.88590404530202224719e+02, -0.92014824378423334394e+02]
cMode[0,3,:] = [-0.95152195796331646704e+02, -0.36956353969242154988e+03, -0.20695307247843501841e+01, -0.38362497103275590148e+03,
-0.22776414558915512031e+03, -0.46624811555373160132e+02, 0.45685596890036768158e+02, -0.25467538299896208542e+03]
cMode[0,4,:] = [0.32918009134216310584e+03, 0.12851583684685703445e+04, 0.71222223072393315845e+01, 0.13249069943127731452e+04,
0.74436590714498462872e+03, 0.14749387222546454623e+03, 0.42968129090712094964e+02, 0.33666974003731491293e+03]
cMode[0,5,:] = [0.20156315877857942098e+03, 0.80514059433962259504e+03, 0.42634222620010087112e+01, 0.79880139299917756190e+03,
0.34272486825284445011e+03, 0.52606372363987690121e+02, 0.74826789775518696146e+01, -0.49206154796489274261e+01]
cMode[0,6,:] = [-0.53421252140832686805e+03, -0.21573846013695305856e+04, -0.11180286753450849879e+02, -0.20919394990520148169e+04,
-0.77083430972020012816e+03, -0.74269446238765679524e+02, 0.86081879160377337001e+01, -0.33189293006747679903e+01]
cMode[0,7,:] = [0.64042667076955872573e+03, 0.25941583461638844099e+04, 0.13364255832056319839e+02, 0.24982768986941823463e+04,
0.87945971890665362025e+03, 0.70070351276938351858e+02, -0.86139878258589916981e+01, 0.30011846979293261838e+01]
cMode[1,0,:] = [-0.63013334128716040893e-13, 0.60865772431920603935e-05, -0.40016262502772708131e-09, -0.23069864565108759713e-04,
0.23255146051530801720e-03,-0.51951144258714805346e-02, 0.36136455088219365805e-01, -0.53017682607232261560e-01]
cMode[1,1,:] = [0.55436122531818750047e-03, -0.66992043381081387565e-02, 0.83803970952918511727e-02, 0.80754210845424374554e-02,
-0.61427751801047527635e-02, 0.54122156026262873140e+00, -0.43010880044939874267e+01, 0.65777678112627633311e+01]
cMode[1,2,:] = [-0.23084376754368056694e+02, 0.20899246010739616075e+03, -0.36836390749521057408e+03, -0.52914996781002994197e+03,
-0.20007855689466294002e+03, -0.36056061195965001253e+02, -0.16132602175333610183e+02, 0.21643026623710106548e+01]
cMode[1,3,:] = [-0.30489331973819799870e+01, 0.33304902492752548326e+02, -0.47908018506183607243e+02, -0.66968634395671120529e+02,
-0.34606055921065168590e+02, -0.74942275618948936966e+01, 0.20818848022147209420e+02, -0.68532587218682348151e+02]
cMode[1,4,:] = [0.44092752992499013586e+01, -0.47413797913825259655e+02, 0.69401964003076273002e+02, 0.97614272148050673649e+02,
0.48704912842393142113e+02, 0.99941378554753921292e+01, -0.14570861805562320689e+02, 0.65068431858518227528e+02]
cMode[1,5,:] = [0.27754841758044097588e+02, -0.24286715291266780525e+03, 0.44379959765035561503e+03, 0.63697225078825798760e+03,
0.23079424266757904149e+03, 0.39146925311929221535e+02, 0.12991883575647105075e+02, -0.34552170740359109402e+01]
cMode[1,6,:] = [-0.20272909693943255149e+02, 0.14757207634388509021e+03, -0.32708628601917539846e+03, -0.46225115463536488036e+03,
-0.13753186496272293837e+03, -0.14358017591372269627e+02, 0.41733787941502891172e+00, -0.14829734668187728452e+00]
cMode[1,7,:] = [0.22054195737479753702e+02, -0.15104824369566767217e+03, 0.35667767962821925742e+03, 0.50036616092734789162e+03,
0.14078440554104824755e+03, 0.12025910667067751802e+02, -0.70328949883151059552e+00, 0.24434145951697048282e+00]
cMode[2,0,:] = [-0.28862183650305684778e-01, 0.16202748718849406373e-02, 0.62452442623701900359e-08, -0.10292552636752272388e-04,
0.14899979598450072693e-10, 0.54701873831866327790e+00, -0.33801096158365382393e+01, 0.46663975792709919687e+01]
cMode[2,1,:] = [0.82086900572783427776e+00, 0.66781995371601139410e+00, 0.17689517609606757453e+01, 0.18458732991639696052e+01,
-0.13408206456871818446e-02, -0.55045820392109447993e+01, 0.34057914707574474810e+02, -0.44172744618314503384e+02]
cMode[2,2,:] = [-0.33440561402073085695e+04, -0.63006668870531044035e+04, -0.16148889754716360123e+05, -0.16734685574071548330e+05,
0.98894432761575998824e+01, -0.10391741179165414621e+04, -0.17857836173455567951e+04, 0.29803489709704852117e+04]
cMode[2,3,:] = [0.31915064954348162373e+04, 0.59270193788343030760e+04, 0.15182397380443129364e+05, 0.15734311019448650625e+05,
-0.93434105903445914265e+01, 0.99260524095183413351e+03, 0.20341329769010583206e+04, -0.32962143597458424260e+04]
cMode[2,4,:] = [-0.92169284592923155230e+02, -0.14814666530740792538e+03, -0.37839662711403616590e+03, -0.39255805350361141492e+03,
0.24567231705288499199e+00, -0.11195431933835682247e+02, -0.27154689478860825069e+03, 0.37730991124653248114e+03]
cMode[2,5,:] = [0.54075863200485683179e+03, 0.13068947201920562140e+04, 0.34015030753731947399e+04, 0.35222027482034166112e+04,
-0.19372427035775936943e+01, 0.10115813738135293053e+03, -0.23832476149210184424e+02, -0.16890808415607296844e+02]
cMode[2,6,:] = [-0.77613909679256050111e+03, -0.22744934351169181496e+04, -0.60171788754328519033e+04, -0.62294563978631467549e+04,
0.32775011290202917813e+01, -0.72131812461142832404e+02, 0.27116603294455088324e+02, -0.71224576605140779150e+01]
cMode[2,7,:] = [0.69369255218888747904e+03, 0.22078310233442990373e+04, 0.58858136079721381506e+04, 0.60932655114670426499e+04,
-0.31519955877359868701e+01, 0.44446156981790601037e+02, -0.14233922505288183479e+02, 0.34387488557126792976e+01]
cMode[3,0,:] = [-0.36868002035455376130e-21, 0.10541748946585167701e-17, 0.27755566547052099579e-13, -0.34780372211895671519e-15,
-0.99328472275559516191e-12, 0.37378357885053432596e-10, -0.27315035342737905565e-09, 0.40965931239852784173e-09]
cMode[3,1,:] = [-0.33737087677076624814e-04, -0.11396625664731094840e-02, -0.10789923154822660400e-01, -0.80796849963401182748e-02,
0.32166227895723737972e-01, -0.25571424523548520468e+01, 0.21267119641439142796e+02, -0.33078925310689939465e+02]
cMode[3,2,:] = [-0.17479871537929662750e+01, -0.72444600080193755076e+02, -0.55033093275560522883e+03, -0.50210145669063086515e+03,
-0.25844521786360195036e+03, -0.55113408951058726614e+02, 0.20187135793924890769e+03, -0.81289982453232934034e+03]
cMode[3,3,:] = [0.78603614447709357904e+00, 0.32329428441630612134e+02, 0.24762820448825531016e+03, 0.22413360241783610860e+03,
0.11996249854833778059e+03, 0.30448969122459486058e+02, -0.23988687682537666034e+03, 0.67671098556149500424e+03]
cMode[3,4,:] = [0.12159934205840530196e+01, 0.51016427583112839982e+02, 0.38232595903016939331e+03, 0.35346430190826185757e+03,
0.17043627550273194870e+03, 0.36343814970370997841e+02, 0.25302217119464871508e+02, 0.16017628158741096910e+03]
cMode[3,5,:] = [-0.13415948717761450037e+01, -0.59032589809690589888e+02, -0.41777999557453684431e+03, -0.40902457469604014406e+03,
-0.14475183395103458749e+03, -0.29478536348031818548e+02, -0.13387094046024579085e+01, 0.23831368497822609242e+01]
cMode[3,6,:] = [0.55621147775238517496e+01, 0.24890779339281028370e+03, 0.17228008107309571883e+04, 0.17256605312262474072e+04,
0.52806415548478540245e+03, 0.64772764588200661961e+02, -0.12923255027738644784e+02, 0.58752559976531326668e+01]
cMode[3,7,:] = [-0.73285857760141963623e+01, -0.32941251542733325230e+03, -0.22662989686584209536e+04, -0.22842752627361360140e+04,
-0.66988125859857134969e+03, -0.67321294156948754405e+02, 0.12094082762406057618e+02, -0.50344380663030463551e+01]
cMode[4,0,:] = [0.57951308146739251014e+01, -0.14999087019121093433e+03, -0.19559326802070165385e+01, -0.66720804039446477418e+03,
-0.75756207610197607849e+02, -0.79695623597154439110e+02, 0.12741669724908954997e+04, -0.25779398744795645193e+04]
cMode[4,1,:] = [0.25559373539070939784e+03, -0.26954324048080380293e+04, 0.42741834334120607508e+02, -0.13609107834915077361e+05,
-0.18219129315199742435e+04, -0.49540503909243787106e+03, -0.86525145704320394202e+03, -0.19349682251463748983e+04]
cMode[4,2,:] = [-0.87923543175638503299e+05, 0.30861003822923382955e+06, -0.25217316373834348652e+05, 0.21333098058144601694e+07,
0.33404584481156045505e+06, 0.42354436002087458845e+05, 0.58637814837743640339e+04, -0.30291614862820930298e+04]
cMode[4,3,:] = [0.21501887419077996277e+06, -0.98416860975998154970e+06, 0.50756266605703403982e+05, -0.58885368624678218196e+07,
-0.88230632458050202160e+06, -0.24671074506661172520e+06, -0.14496629047056444505e+06, 0.68193611837692298394e+05]
cMode[4,4,:] = [0.18016800297232990146e+02, 0.11854731605905291047e+00, 0.41242056558479642802e-02, 0.10788775915866732901e+00,
0.76539029541206167195e+00, 0.67292188204585867694e+02, -0.31151519702859546967e+03, 0.44202304203357769551e+03]
cMode[4,5,:] = [0.20068140169309098830e+06, -0.93099363016509286694e+06, 0.47215756436253171202e+05, -0.55473705494319487385e+07,
-0.82963829864856943885e+06, -0.23595446376207669381e+06, -0.14387731979865256981e+06, 0.66726697732453104094e+05]
cMode[4,6,:] = [-0.10659833147754620430e+06, 0.36525155976544767533e+06, -0.32972599682804299980e+05, 0.26280915665374022793e+07,
0.41399542203693986408e+06, 0.40200322950430580348e+05, 0.33714384362519780324e+04, -0.16170695162643912823e+04]
cMode[4,7,:] = [0.45385852289903034773e+05, -0.15367285853252436567e+06, 0.15938904084451054998e+05, -0.11759967065602889846e+07,
-0.18617563074881875451e+06, -0.11629254312624910383e+05, -0.37356632983270783299e+03, 0.13368622059973152005e+03]
cMode[5,0,:] = [0.37071935082976326114e-15, 0.46201169190899653571e-08, -0.62075147292222574435e-11, -0.36146708876274598054e-06,
0.87046861427880060091e-05, -0.10076975633470059978e-03, 0.89298240723925399464e-03, -0.14060364310276510124e-02]
cMode[5,1,:] = [-0.56188781974792609830e-01, -0.40650854381275376425e+01, -0.14988738284605299000e+01, -0.38632728626301734209e+01,
-0.22261011246466999580e+01, -0.14212548821848476343e+02, 0.19383430976405556123e+03, -0.34806064213022160913e+03]
cMode[5,2,:] = [0.27593841470781566016e+03, 0.22140385471263797079e+05, 0.78362796461845718454e+04, 0.19668830633326862766e+05,
0.10181292569179314355e+05, 0.39096643822536609746e+04, 0.26876686999693699675e+04, 0.33310618692009836827e+04]
cMode[5,3,:] = [0.44061551736067183782e+02, 0.34746817101944778016e+04, 0.12384774155446498511e+04, 0.31312056838986070950e+04,
0.17260274467090328087e+04, 0.70922681375271823256e+03, -0.30968222241454737009e+03, 0.29976750308301927105e+04]
cMode[5,4,:] = [-0.29826526924997253331e+03, -0.23828019858653441964e+05, -0.84485266842828430355e+04, -0.21247202725657969857e+05,
-0.11181459854792739072e+05, -0.43708622712852429614e+04, -0.24560179552632419586e+04, -0.60645360750106398484e+04]
cMode[5,5,:] = [-0.74893303991611572811e+02, -0.62909773314430736945e+04, -0.21848887114479706994e+04, -0.53368828632132494504e+04,
-0.22162146435679854761e+04, -0.52579256909996399116e+03, -0.11309723650631768876e+03, 0.63863048387678587047e+02]
cMode[5,6,:] = [0.19717911116782794067e+03, 0.16899744672548959734e+05, 0.58199397638247170050e+04, 0.13983841291452328015e+05,
0.50676093636595025415e+04, 0.75432123378980122652e+03, 0.10996364977646411276e+02, -0.48990730626876937137e+01]
cMode[5,7,:] = [-0.22796489432205970793e+03, -0.19665709629971770411e+05, -0.67538957248075544015e+04, -0.16130393271007925193e+05,
-0.55557948391593710013e+04, -0.67829337377268315023e+03, 0.63691256079676055179e+01, -0.36820923261738309761e+01]
####
#### ; coeff p
#### p[0,*] = [0.32901733177624610249d-01, 0.54475681390900216882d-01, 0.11033542666255162778d+00, $
#### 0.21824017019698205288d+00, 0.43540342906595910221d+00, 0.87057139406808232706d+00, $
#### 0.17411041057694376377d+01, 0.34822027229232799250d+01, 0.29391078654639329670d-01, $
#### 0.53092570701901200536d-01, 0.10068022495463000431d+00, 0.21707901052684266396d+00, $
#### 0.43526892089758080217d+00, 0.87054714666225887498d+00, 0.17411005963420354447d+01, $
#### 0.34822023079824373503d+01]
#### p[1,*] = [0.32535198971349430507d-01, 0.57792553225847100861d-01, 0.11708104989648597804d+00, $
#### 0.21135795367821805790d+00, 0.42489217791799518408d+00, 0.86257818103884389415d+00, $
#### 0.17181130605554559842d+01, 0.34339371403040641617d+01, 0.25660275820190880935d-01, $
#### 0.45569108067616728163d-01, 0.86363227093303382986d-01, 0.20636215029847653212d+00, $
#### 0.42568882948609800820d+00, 0.85912749746654117899d+00, 0.17163445496343952001d+01, $
#### 0.34350464850833977159d+01]
#### p[2,*] = [0.32114531391568896800d-01, 0.54668167014354631660d-01, 0.10999415084820487464d+00, $
#### 0.21830735865316683863d+00, 0.43522952245694588313d+00, 0.87059254539023882557d+00, $
#### 0.17410412227916809868d+01, 0.34821651743153156921d+01, 0.25685870467860096866d-01, $
#### 0.45698157579046663201d-01, 0.90356977415432240263d-01, 0.21884901172766766386d+00, $
#### 0.43511983572804924236d+00, 0.87029548334738358050d+00, 0.17410538620391950992d+01, $
#### 0.34822549594214162738d+01]
#### p[3,*] = [0.32356337372804593321d-01, 0.54152800183054496940d-01, 0.10961184933750356407d+00, $
#### 0.21794853515003347332d+00, 0.43533813536787366871d+00, 0.87056003928681793269d+00, $
#### 0.17411018194244403112d+01, 0.34822024269907854154d+01, 0.27588798509226557520d-01, $
#### 0.48812593148800704767d-01, 0.93213066985517798457d-01, 0.21740973446223286202d+00, $
#### 0.43526939652742839825d+00, 0.87054487704825334049d+00, 0.17411007942523434977d+01, $
#### 0.34822023168779896451d+01]
#### p[4,*] = [0.32735344319694603676d-01, 0.54347222488222799441d-01, 0.10989084818706484902d+00, $
#### 0.21805120636819852464d+00, 0.43536325140487690532d+00, 0.87056578462259590622d+00, $
#### 0.17411028067887455605d+01, 0.34822028139735676788d+01, 0.27779830466381087994d-01, $
#### 0.50806515286805851161d-01, 0.95276608648916276678d-01, 0.21745769582596392588d+00, $
#### 0.43527860088967376128d+00, 0.87054808117661632849d+00, 0.17411000172672461694d+01, $
#### 0.34822021788641484008d+01]
#### p[5,*] = [0.33568077258163553366d-01, 0.53905902840230535133d-01, 0.10958409128391679487d+00, $
#### 0.21675427898182215713d+00, 0.43240630217346920360d+00, 0.86456015794363612059d+00, $
#### 0.17290777196751532684d+01, 0.34581494323600296958d+01, 0.28810823097181398111d-01, $
#### 0.53751387895676616679d-01, 0.10088751776373805491d+00, 0.21567365023402529367d+00, $
#### 0.43226560363373209838d+00, 0.86453474685580449232d+00, 0.17290747052830268692d+01, $
#### 0.34581487532588752742d+01]
p[0,:] = [0.32901733177624610249e-01, 0.54475681390900216882e-01, 0.11033542666255162778e+00,
0.21824017019698205288e+00, 0.43540342906595910221e+00, 0.87057139406808232706e+00,
0.17411041057694376377e+01, 0.34822027229232799250e+01, 0.29391078654639329670e-01,
0.53092570701901200536e-01, 0.10068022495463000431e+00, 0.21707901052684266396e+00,
0.43526892089758080217e+00, 0.87054714666225887498e+00, 0.17411005963420354447e+01,
0.34822023079824373503e+01]
p[1,:] = [0.32535198971349430507e-01, 0.57792553225847100861e-01, 0.11708104989648597804e+00,
0.21135795367821805790e+00, 0.42489217791799518408e+00, 0.86257818103884389415e+00,
0.17181130605554559842e+01, 0.34339371403040641617e+01, 0.25660275820190880935e-01,
0.45569108067616728163e-01, 0.86363227093303382986e-01, 0.20636215029847653212e+00,
0.42568882948609800820e+00, 0.85912749746654117899e+00, 0.17163445496343952001e+01,
0.34350464850833977159e+01]
p[2,:] = [0.32114531391568896800e-01, 0.54668167014354631660e-01, 0.10999415084820487464e+00,
0.21830735865316683863e+00, 0.43522952245694588313e+00, 0.87059254539023882557e+00,
0.17410412227916809868e+01, 0.34821651743153156921e+01, 0.25685870467860096866e-01,
0.45698157579046663201e-01, 0.90356977415432240263e-01, 0.21884901172766766386e+00,
0.43511983572804924236e+00, 0.87029548334738358050e+00, 0.17410538620391950992e+01,
0.34822549594214162738e+01]
p[3,:] = [0.32356337372804593321e-01, 0.54152800183054496940e-01, 0.10961184933750356407e+00,
0.21794853515003347332e+00, 0.43533813536787366871e+00, 0.87056003928681793269e+00,
0.17411018194244403112e+01, 0.34822024269907854154e+01, 0.27588798509226557520e-01,
0.48812593148800704767e-01, 0.93213066985517798457e-01, 0.21740973446223286202e+00,
0.43526939652742839825e+00, 0.87054487704825334049e+00, 0.17411007942523434977e+01,
0.34822023168779896451e+01]
p[4,:] = [0.32735344319694603676e-01, 0.54347222488222799441e-01, 0.10989084818706484902e+00,
0.21805120636819852464e+00, 0.43536325140487690532e+00, 0.87056578462259590622e+00,
0.17411028067887455605e+01, 0.34822028139735676788e+01, 0.27779830466381087994e-01,
0.50806515286805851161e-01, 0.95276608648916276678e-01, 0.21745769582596392588e+00,
0.43527860088967376128e+00, 0.87054808117661632849e+00, 0.17411000172672461694e+01,
0.34822021788641484008e+01]
p[5,:] = [0.33568077258163553366e-01, 0.53905902840230535133e-01, 0.10958409128391679487e+00,
0.21675427898182215713e+00, 0.43240630217346920360e+00, 0.86456015794363612059e+00,
0.17290777196751532684e+01, 0.34581494323600296958e+01, 0.28810823097181398111e-01,
0.53751387895676616679e-01, 0.10088751776373805491e+00, 0.21567365023402529367e+00,
0.43226560363373209838e+00, 0.86453474685580449232e+00, 0.17290747052830268692e+01,
0.34581487532588752742e+01]
####
#### ; coeff r
#### r[0,*] = [0.23400582376918621640d-01, 0.29906910603763133593d+00, 0.59419217080359816307d+00, $
#### 0.51683309796288279258d+00, 0.55280210585007472090d+00, 0.89024092184836405294d+00, $
#### 0.17236578859735887547d+01, 0.34342717240262339295d+01, 0.89287729231317083389d-01, $
#### 0.15439963383324775136d+00, 0.50639243810157408276d-01, 0.24456407532792816539d+00, $
#### 0.39627924482769603997d+00, 0.86233612861692066076d+00, 0.17188649618875711411d+01, $
#### 0.34339261899212027984d+01]
#### r[1,*] = [0.70036641254021372304d-01, 0.33297901220247387854d+00, 0.80631162662909527938d+00, $
#### 0.54483599484759341891d+00, 0.56290677049809243470d+00, 0.87422019588133075274d+00, $
#### 0.17295204193052738261d+01, 0.34651766994561397083d+01, 0.57115996571135090320d-01, $
#### 0.23559291883706992010d+00, 0.11690775021980281955d+00, 0.26929992715503474620d+00, $
#### 0.42836762405666322095d+00, 0.87049839288489376798d+00, 0.17350924596759142559d+01, $
#### 0.34557476257228634253d+01]
#### r[2,*] = [0.11127103098086079668d+00, 0.42445608810651771491d+00, 0.71058432225687280237d+00, $
#### 0.69760620618057709307d+00, 0.58851822769682966551d+00, 0.10410712531090018373d+01, $
#### 0.18016710591015945297d+01, 0.35162043268667866335d+01, 0.50015594405195891170d+00, $
#### 0.32951760674056020938d+00, 0.11251791947458271714d+00, 0.19597524022502758711d+00, $
#### -0.52856178098938624287d-01, 0.95266905942729742662d+00, 0.18279298809137124237d+01, $
#### 0.35144044692591625000d+01]
#### r[3,*] = [0.22860180522694295568d-01, 0.25001972739726898709d+00, 0.46402414947596950511d+00, $
#### 0.44398927627952060603d+00, 0.50272011857581864191d+00, 0.88445628872250843244d+00, $
#### 0.17389166922733014786d+01, 0.34834653146383227628d+01, 0.49049989399817732760d-01, $
#### 0.85147557568237122183d-01, 0.23424821227315173466d+00, 0.14250169643289223309d+00, $
#### 0.40586458823694400166d+00, 0.87361464455904407344d+00, 0.17279679104606691097d+01, $
#### 0.34828619038557646625d+01]
#### r[4,*] = [0.48558255088791124620d+00, -0.60422350819111194653d+00,-0.13683535645881987896d+01, $
#### -0.85855064913451286656d+00, 0.28100294887963377377d+00, 0.84995521582071482669d+00, $
#### 0.17002422172314077819d+01, 0.34601959118742100507d+01, -0.60269701842058500673d+00, $
#### -0.31233918670609019940d+00, -0.71139286069860601102d-01, 0.13800610282920311533d+00, $
#### 0.38218789451073229557d+00, 0.81980250822178710734d+00, 0.17183601833584307705d+01, $
#### 0.34710558209881763325d+01]
#### r[5,*] = [0.38821262813796547419d-01, 0.35983028187700409894d+00, 0.57018233005798357737d+00, $
#### 0.51153414239116576922d+00, 0.55367851931772325002d+00, 0.92841885168264770555d+00, $
#### 0.17036409347489702703d+01, 0.34382652552776309384d+01, 0.54101616716924096905d-01, $
#### 0.17168542829213128797d+00, 0.10028211367824637623d+00, 0.27153034480893150082d+00, $
#### 0.43311645903980116045d+00, 0.75314774033930556029d+00, 0.17038798689307688150d+01, $
#### 0.34370447884477832722d+01]
r[0,:] = [0.23400582376918621640e-01, 0.29906910603763133593e+00, 0.59419217080359816307e+00,
0.51683309796288279258e+00, 0.55280210585007472090e+00, 0.89024092184836405294e+00,
0.17236578859735887547e+01, 0.34342717240262339295e+01, 0.89287729231317083389e-01,
0.15439963383324775136e+00, 0.50639243810157408276e-01, 0.24456407532792816539e+00,
0.39627924482769603997e+00, 0.86233612861692066076e+00, 0.17188649618875711411e+01,
0.34339261899212027984e+01]
r[1,:] = [0.70036641254021372304e-01, 0.33297901220247387854e+00, 0.80631162662909527938e+00,
0.54483599484759341891e+00, 0.56290677049809243470e+00, 0.87422019588133075274e+00,
0.17295204193052738261e+01, 0.34651766994561397083e+01, 0.57115996571135090320e-01,
0.23559291883706992010e+00, 0.11690775021980281955e+00, 0.26929992715503474620e+00,
0.42836762405666322095e+00, 0.87049839288489376798e+00, 0.17350924596759142559e+01,
0.34557476257228634253e+01]
r[2,:] = [0.11127103098086079668e+00, 0.42445608810651771491e+00, 0.71058432225687280237e+00,
0.69760620618057709307e+00, 0.58851822769682966551e+00, 0.10410712531090018373e+01,
0.18016710591015945297e+01, 0.35162043268667866335e+01, 0.50015594405195891170e+00,
0.32951760674056020938e+00, 0.11251791947458271714e+00, 0.19597524022502758711e+00,
-0.52856178098938624287e-01, 0.95266905942729742662e+00, 0.18279298809137124237e+01,
0.35144044692591625000e+01]
r[3,:] = [0.22860180522694295568e-01, 0.25001972739726898709e+00, 0.46402414947596950511e+00,
0.44398927627952060603e+00, 0.50272011857581864191e+00, 0.88445628872250843244e+00,
0.17389166922733014786e+01, 0.34834653146383227628e+01, 0.49049989399817732760e-01,
0.85147557568237122183e-01, 0.23424821227315173466e+00, 0.14250169643289223309e+00,
0.40586458823694400166e+00, 0.87361464455904407344e+00, 0.17279679104606691097e+01,
0.34828619038557646625e+01]
r[4,:] = [0.48558255088791124620e+00, -0.60422350819111194653e+00,-0.13683535645881987896e+01,
-0.85855064913451286656e+00, 0.28100294887963377377e+00, 0.84995521582071482669e+00,
0.17002422172314077819e+01, 0.34601959118742100507e+01, -0.60269701842058500673e+00,
-0.31233918670609019940e+00, -0.71139286069860601102e-01, 0.13800610282920311533e+00,
0.38218789451073229557e+00, 0.81980250822178710734e+00, 0.17183601833584307705e+01,
0.34710558209881763325e+01]
r[5,:] = [0.38821262813796547419e-01, 0.35983028187700409894e+00, 0.57018233005798357737e+00,
0.51153414239116576922e+00, 0.55367851931772325002e+00, 0.92841885168264770555e+00,
0.17036409347489702703e+01, 0.34382652552776309384e+01, 0.54101616716924096905e-01,
0.17168542829213128797e+00, 0.10028211367824637623e+00, 0.27153034480893150082e+00,
0.43311645903980116045e+00, 0.75314774033930556029e+00, 0.17038798689307688150e+01,
0.34370447884477832722e+01]
#### pvec = (rvec = dblarr(16))
pvec = np.array([16])
rvec = np.array([16])
#### ain = (cin = dblarr(M,M))
ain = np.array((M,M))
cin = np.array((M,M))
#### CS = [1.068627d, 1.572350d, 1.015095d, 1.714487d, -5.232800d, 6.528898d]/0.002d
CS = np.array([1.068627, 1.572350, 1.015095, 1.714487, -5.232800, 6.528898])/0.002
####
#### Bx = 0d
#### By = 0d
#### Bz = 0d
Bx = 0.
By = 0.
Bz = 0.
####
#### Bxout = 0d
#### Byout = 0d
#### Bzout = 0d
Bxout = 0.
Byout = 0.
Bzout = 0.
####
#### xp = [x + 0.001d, x, x]
#### yp = [y, y + 0.001d, y]
#### zp = [z, z, z + 0.001d]
xp = [x + 0.001, x, x]
yp = [y, y + 0.001, y]
zp = [z, z, z + 0.001]
#### xm = [x - 0.001d, x, x]
#### ym = [y, y - 0.001d, y]
#### zm = [z, z, z - 0.001d]
xm = [x - 0.001, x, x]
ym = [y, y - 0.001, y]
zm = [z, z, z - 0.001]
#### for Mode = 0, 5 do begin
for Mode in range(6):
#### ain = reform(aMode[Mode,*,*])
ain = aMode[Mode,:,:]
#### cin = reform(cMode[Mode,*,*])
cin = cMode[Mode,:,:]
#### p_1 = reform(p[Mode,0:7])
p_1 = p[Mode,0:8]
#### p_2 = reform(p[Mode,8:15])
p_2 = p[Mode,8:16]
#### r_1 = reform(r[Mode,0:7])
r_1 = r[Mode,0:8]
#### r_2 = reform(r[Mode,8:15])
r_2 = r[Mode,8:16]
####
#### Up = dblarr(3)
#### Um = Up
Up = np.zeros([3])
Um = Up
#### term_p = (term_r = dblarr(8,8))
term_p = np.zeros([8, 8])
term_r = np.zeros([8, 8])
#### for i = 0, 7 do begin
#### term_p[i,0:7] = 1d/(p_1[i]*p_1[i]) + 1d/(p_2[0:7]*p_2[0:7])
#### term_r[i,0:7] = 1d/(r_1[i]*r_1[i]) + 1d/(r_2[0:7]*r_2[0:7])
#### endfor
for i in range(8):
term_p[i,0:8] = 1./(p_1[i]*p_1[i]) + 1./(p_2[0:8]*p_2[0:8])
term_r[i,0:8] = 1./(r_1[i]*r_1[i]) + 1./(r_2[0:8]*r_2[0:8])
#### term_p = exp(sqrt(term_p))
#### term_r = exp(sqrt(term_r))
term_p = np.exp(np.sqrt(term_p))
term_r = np.exp(np.sqrt(term_r))
####
#### for i = 0, 7 do begin
for i in range(8):
#### cos_yp = cos(yp/p_1[i])
cos_yp = np.cos(yp/p_1[i])
#### cos_ym = cos(ym/p_1[i])
cos_ym = np.cos(ym/p_1[i])
#### sin_yp = sin(yp/r_1[i])
sin_yp = np.sin(yp/r_1[i])
#### sin_ym = sin(ym/r_1[i])
sin_ym = np.sin(ym/r_1[i])
#### for k = 0, 7 do begin
for k in range(8):
#### Up = Up + ain[i,k]*(term_p[i,k]^xp)*cos_yp*sin(zp/p_2[k]) + cin[i,k]*(term_r[i,k]^xp)*sin_yp*sin(zp/r_2[k])
#### Um = Um + ain[i,k]*(term_p[i,k]^xm)*cos_ym*sin(zm/p_2[k]) + cin[i,k]*(term_r[i,k]^xm)*sin_ym*sin(zm/r_2[k])
Up = Up + ain[i,k]*(term_p[i,k] ** xp)*cos_yp*np.sin(zp/p_2[k]) + cin[i,k]*(term_r[i,k] ** xp)*sin_yp*np.sin(zp/r_2[k])
Um = Um + ain[i,k]*(term_p[i,k] ** xm)*cos_ym*np.sin(zm/p_2[k]) + cin[i,k]*(term_r[i,k] ** xm)*sin_ym*np.sin(zm/r_2[k])
#### endfor
#### endfor
#### B = Up - Um
B = Up - Um
#### Bxout += B[0]*CS[Mode]
#### Byout += B[1]*CS[Mode]
#### Bzout += B[2]*CS[Mode]
Bxout += B[0]*CS[Mode]
Byout += B[1]*CS[Mode]
Bzout += B[2]*CS[Mode]
#### endfor
### idl
# -0.17749937 -0.023122741 -2.8439913
### Python
# DEBUG:kk2009:Bxout=-0.17749938914011432
# DEBUG:kk2009:Byout=-0.023122815271786192
# DEBUG:kk2009:Bzout=-2.843991083303587
### Sofarsogood, 2014-05-04T19:11:50
_logger.debug('Bxout={}'.format(Bxout))
_logger.debug('Byout={}'.format(Byout))
_logger.debug('Bzout={}'.format(Bzout))
return Bxout, Byout, Bzout
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure tail_mag_shield_cyl_S3
#### ;--------------------------------------------------------------------------------------
#### pro tail_mag_shield_cyl_S3, ctime, rmap, pmap, zmap, Brcss, Bpcss, Bzcss, M, sphiOut, Mode
[docs]def tail_mag_shield_cyl_S3(ctime, rmap, pmap, zmap, M, sphiOut, Mode):
#### ; transform to cartesian
#### cyl2car_pos, rmap, pmap, zmap, Xcar, Ycar, Zcar
Xcar, Ycar, Zcar = cyl2car_pos(rmap, pmap, zmap)
####
#### S32PJSO, xJSO, yJSO, zJSO, Xcar, Ycar, Zcar, sphiOut
xJSO, yJSO, zJSO = S32PJSO(Xcar, Ycar, Zcar, sphiOut)
####
#### ; Divide position by 100 to match the coefficients
#### ; checkIfInsideMagnetop, xJSO, yJSO, zJSO, ianswer
#### xJSO = xJSO/100d
#### yJSO = yJSO/100d
#### zJSO = zJSO/100d
xJSO = xJSO/100.
yJSO = yJSO/100.
zJSO = zJSO/100.
####
#### ; works in the original coordinate system S3 ***
#### B_tail_shield, xJSO, yJSO, zJSO, BxJSO, ByJSO, BzJSO, M, Mode
BxJSO, ByJSO, BzJSO = B_tail_shield(xJSO, yJSO, zJSO, M, Mode)
####
#### PJSO2S3, BxJSO, ByJSO, BzJSO, BxS3, ByS3, BzS3, sphiOut
BxS3, ByS3, BzS3 = PJSO2S3(BxJSO, ByJSO, BzJSO, sphiOut)
####
#### ; transform back to cylindrical
#### car2cyl_vect, Brcss, Bpcss, Bzcss, BxS3, ByS3, BzS3, pmap
Brcss, Bpcss, Bzcss = car2cyl_vect(BxS3, ByS3, BzS3, pmap)
return Brcss, Bpcss, Bzcss
#### end ; tail_mag_shield_cyl_S3
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure csheet_struc
#### ;
#### ; NEW VERSION NOV 2003
#### ; This program calculates the distance of current sheet from the
#### ; Jovigraphic equator in system III coordinates
#### ;--------------------------------------------------------------------------------------
#### pro csheet_struc, ZNS3, rho, phi, XJSO, YJSO, stheta, ctime ; ctime?
[docs]def csheet_struc(rho, phi, XJSO, YJSO, stheta, ctime):
''' This program calculates the distance of current sheet from the Jovigraphic equator in system III coordinates
>>> zns3 = csheet_struc(30.0, 0.5, 28.0, 1.2, np.deg2rad(20.), ctimer(2000, 1, 1, 0, 0, 0))
>>> print('%.5f' % zns3)
2.54079
The returned value is examined by IDL.
::
IDL> .compile kk_2009
IDL> csheet_struc, retval, 30.0, 0.5, 28.0, 1.2, 20. * 3.1416 / 180., ctimer(2000, 1, 1, 0, 0, 0)
IDL> print, retval
2.5407947
'''
#### radian = 0.01745329252d
radian = 0.01745329252
#### period = 9.927953d
period = 9.927953
#### OMEGAJ = 36.26125d
OMEGAJ = 36.26125
#### X0 = -45d
X0 = -45.
#### phip0 = 6.12611d
phip0 = 6.12611
#### ;ten = tan(0.167551606d)
#### ten = 0.16913733749783d
ten = 0.16913733749783
#### C1 = 0.005973d
C1 = 0.005973
#### C2 = 5.114d-5
C2 = 5.114e-5
#### c3 = 1.59d-5
C3 = 1.59e-5
#### psi2 = -1.201201d
psi2 = -1.201201
#### C4 = 0.313244d
C4 = 0.313244
#### C5 = -0.366166d
C5 = -0.366166
#### psi4 = 2.2604522d
psi4 = 2.2604522
####
#### if (abs(XJSO) lt 1d-6) then XJSO = 1d-6
if abs(XJSO) < 1e-6:
XJSO = 1e-6 #YF: Probably to avoid singular point...
####
#### RLT = atan(YJSO, XJSO) ; calculate RLT here
RLT = np.arctan2(YJSO, XJSO) # calculate RLT here
####
#### delay1 = C1*rho + 0.5d*C2*(rho*rho)*cos(RLT - psi2) + 0.5d*C3*rho*rho
delay1 = C1*rho + 0.5*C2*(rho*rho)*np.cos(RLT - psi2) + 0.5*C3*rho*rho
#### Alfven = -360d/( period*(C4*cos(RLT - psi4) + C5) )
Alfven = -360/( period*(C4*np.cos(RLT - psi4) + C5) )
#### delay2 = rho/Alfven*Omegaj*radian
delay2 = rho/Alfven*OMEGAJ*radian
####
#### phip = phip0 - delay1 - delay2
phip = phip0 - delay1 - delay2
####
#### hyptan = x0*tanh(XJSO/x0)
hyptan = X0*np.tanh(XJSO/X0)
#### rho1 = sqrt(hyptan*hyptan + YJSO*YJSO)
rho1 = np.sqrt(hyptan*hyptan + YJSO*YJSO)
####
#### ZNS3 = rho1*ten*cos(phi - phip) + XJSO*(1d - tanh( abs(x0/xjso) ))*tan(stheta)
ZNS3 = rho1*ten*np.cos(phi - phip) + XJSO*(1. - np.tanh( abs(X0/XJSO) ))*np.tan(stheta)
return ZNS3
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure csheet_N
#### ;
#### ; This procedure uses csheet_structure to calculate the normal
#### ; direction to the current sheet in system III coordinates
#### ;--------------------------------------------------------------------------------------
############# YF considers that RNi is returned, and the other input
[docs]def csheet_N(rho, phi, zs3, stheta, ctime):
r''' Something related to current sheet...
:param rho: Distance from Jupiter (?)
:param zs3: z-value in system III.
:param phi: Phase angle (longitude) in system III, radians.
:param stheta: Angle (radian) in system III(?)
:param ctime: Time in J2000 epoch (?), where anyway :func:`ctimer` derives.
.. note::
In the source code,
.. math::
xs3 = \rho * \cos(\phi)
is done. So that ``phi`` is the angle in the same coordinate system to ``xs3`` (and likely ``zs3``).
Later, the ``rho``, ``phi`` based vector (``xs3``, ``ys3``, ``zs3``) is converted from System3
to JSO. Therefore, they looks in SYS3 coordinate.
>>> ctime = ctimer(2000, 6, 1, 10, 32, 17.0)
>>> print(ctime)
1086085937.0
>>> rnx, rny, rnz = csheet_N(29.955566, 2.7565556, 1.6321992, 0.05433560, ctime)
>>> print('{:.6f} {:.6f} {:.6f}'.format(rnx, rny, rnz))
-0.051022 0.055858 0.997134
.. note::
In IDL original program, the results will be
(-0.050987206 0.055844548 0.99713675)
from the same input.
The difference (only 3-places accuracy) could come from the
fact that the calculation in IDL was done in single-precision,
where 7-places accuracy is there, but the part of numerical
differentiation, similar values are subtracted and only
3-places accuray is left (cancellation of significant digits).
'''
#### pro csheet_N, RNx, RNy, RNz, rho, phi, zs3, stheta, ctime
#### vecin = (vecou = dblarr(3))
vecin = np.zeros([3])
vecou = np.zeros([3])
#### delta = 0.1d
delta = 0.1 ### it looks this function calculate "numerical differentials"
#### delta2 = 2d*delta
delta2 = delta * 2.
#### xs3 = rho*cos(phi)
xs3 = rho * np.cos(phi)
#### ys3 = rho*sin(phi)
ys3 = rho * np.sin(phi)
#### xp = xs3 + delta
xp = xs3 + delta
#### xm = xs3 - delta
xm = xs3 - delta
#### yp = ys3 + delta
yp = ys3 + delta
#### ym = ys3 - delta
ym = ys3 - delta
####
#### ; First calculate the x derivatives
#### ; Define the positive parameters
#### rhop = sqrt(xp*xp + ys3*ys3)
rhop = np.sqrt(xp ** 2 + ys3 ** 2)
#### phip = atan(ys3, xp)
phip = np.arctan2(ys3, xp)
#### vecin = [xp, ys3, zs3]
vecin = np.array([xp, ys3, zs3])
#### JROT, 's3c', 'jso', vecin, vecou, ctime
vecou = jrot('s3c', 'jso', vecin, ctime)
#### XJSOP = vecou[0]
#### YJSOP = vecou[1]
xjsop, yjsop = vecou[:2]
#### csheet_struc, ZNS3P, rhop, phip, XJSOP, YJSOP, stheta, ctime
zns3p = csheet_struc(rhop, phip, xjsop, yjsop, stheta, ctime)
_logger.debug('zns3p = {}'.format(zns3p))
#### ; Now define the negative parameters
#### rhom = sqrt(xm*xm + ys3*ys3)
rhom = np.sqrt(xm*xm + ys3*ys3)
#### phim = atan(ys3,xm)
phim = np.arctan2(ys3,xm)
####
#### vecin = [xm, ys3, zs3]
vecin = np.array([xm, ys3, zs3])
#### JROT, 's3c', 'jso', vecin, vecou, ctime
vecou = jrot('s3c', 'jso', vecin, ctime)
#### XJSOM = vecou[0]
#### YJSOM = vecou[1]
xjsom, yjsom = vecou[:2]
#### csheet_struc, ZNS3M, rhom, phim, XJSOM, YJSOM, stheta, ctime
zns3m = csheet_struc(rhom, phim, xjsom, yjsom, stheta, ctime)
_logger.debug('zns3m = {}'.format(zns3m))
#### dzdx = (ZNS3P - ZNS3M)/delta2
dzdx = (zns3p - zns3m)/delta2
_logger.debug('dzdx = {}'.format(dzdx))
####
#### ; Next the y derivative
#### ; Define the positive parameters
#### rhop = sqrt(xs3*xs3 + yp*yp)
#### phip = atan(yp, xs3)
rhop = np.sqrt(xs3*xs3 + yp*yp)
phip = np.arctan2(yp, xs3)
#### vecin = [xs3, yp, zs3]
vecin = [xs3, yp, zs3]
#### JROT, 's3c', 'jso', vecin, vecou, ctime
vecou = jrot('s3c', 'jso', vecin, ctime)
#### XJSOP = vecou[0]
#### YJSOP = vecou[1]
xjsop = vecou[0]
yjsop = vecou[1]
#### csheet_struc, ZNS3P, rhop, phip, XJSOP, YJSOP, stheta, ctime
zns3p = csheet_struc(rhop, phip, xjsop, yjsop, stheta, ctime)
_logger.debug('zns3p = {}'.format(zns3p))
#### ; Now define the negative parameters
#### rhom = sqrt(xs3*xs3 + ym*ym)
#### phim = atan(ym, xs3)
rhom = np.sqrt(xs3*xs3 + ym*ym)
phim = np.arctan2(ym, xs3)
#### vecin = [xs3, ym, zs3]
vecin = [xs3, ym, zs3]
#### JROT, 's3c', 'jso', vecin, vecou, ctime
vecou = jrot('s3c', 'jso', vecin, ctime)
#### XJSOM = vecou[0]
#### YJSOM = vecou[1]
xjsom = vecou[0]
yjsom = vecou[1]
####
#### csheet_struc, ZNS3M, rhom, phim, XJSOM, YJSOM, stheta, ctime
zns3m = csheet_struc(rhom, phim, xjsom, yjsom, stheta, ctime)
_logger.debug('zns3m = {}'.format(zns3m))
#### dzdy = (ZNS3P - ZNS3M)/delta2
dzdy = (zns3p - zns3m)/delta2
_logger.debug('dzdy = {}'.format(dzdy))
#### RNx = -dzdx
#### RNy = -dzdy
#### RNz = 1d
RNx = -dzdx
RNy = -dzdy
RNz = 1.
#### RN = sqrt(RNx*RNx + RNy*RNy + RNz*RNz)
RN = np.sqrt(RNx*RNx + RNy*RNy + RNz*RNz)
RNx = RNx/RN
RNy = RNy/RN
RNz = RNz/RN
return np.array([RNx, RNy, RNz])
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure get_mapped_sunangle
#### ;--------------------------------------------------------------------------------------
#### pro get_mapped_sunangle, sthetaIN, sphiIN, sthetaOut, sphiOut, $
#### xpx, xpy, xpz, ypx, ypy, ypz, zpx, zpy, zpz
[docs]def get_mapped_sunangle(sthetaIN, sphiIN, xpx, xpy, xpz, ypx, ypy, ypz, zpx, zpy, zpz):
####
#### SPH2CAR_pos, Xin, Yin, Zin, 1d, sthetaIN, sphiIN
Xin, Yin, Zin = SPH2CAR_pos(1., sthetaIN, sphiIN)
####
#### Xout = Xin*xpx + Yin*xpy + Zin*xpz
#### Yout = Xin*ypx + Yin*ypy + Zin*ypz
#### Zout = Xin*zpx + Yin*zpy + Zin*zpz
Xout = Xin*xpx + Yin*xpy + Zin*xpz
Yout = Xin*ypx + Yin*ypy + Zin*ypz
Zout = Xin*zpx + Yin*zpy + Zin*zpz
####
#### CAR2SPH_pos, Xout, Yout, Zout, Rnotused, sthetaOut, sphiOut
Rnotused, sthetaOut, sphiOut = CAR2SPH_pos(Xout, Yout, Zout)
return sthetaOut, sphiOut
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure getBimf
#### ;
#### ; Modified by Fran Bagenal Nov. 2008 - input background field
#### ; Br=Bx = +-0.2, Bp=By = -+1.0, Bt=Bz = 0.0 - in nT - OR Zero...
#### ;--------------------------------------------------------------------------------------
#### pro getBimf, ctime, thetaS3, phiS3, BrS3IMF, BtS3IMF, BpS3IMF
[docs]def getBimf(ctime, thetaS3, phiS3):
#### vecin = (vecou = dblarr(3))
vecin = np.zeros([3])
vecou = np.zeros([3])
#### BxIMFgsm = 0d
#### ByIMFgsm = 0d
#### BzIMFgsm = 0d
BxIMFgsm = 0.
ByIMFgsm = 0.
BzIMFgsm = 0.
####
#### ; go from JSM to S3
#### vecin = [BxIMFgsm, ByIMFgsm, BzIMFgsm]
vecin = [BxIMFgsm, ByIMFgsm, BzIMFgsm]
#### JROT, "jsm", "s3c", vecin, vecou, ctime
vecou = jrot("jsm", "s3c", vecin, ctime)
#### BxIMFs3 = vecou[0]
#### ByIMFs3 = vecou[1]
#### BzIMFs3 = vecou[2]
BxIMFs3 = vecou[0]
ByIMFs3 = vecou[1]
BzIMFs3 = vecou[2]
####
#### ; go to spherical
BrS3IMF, BtS3IMF, BpS3IMF, = CAR2SPH_MAG(BxIMFs3, ByIMFs3, BzIMFs3, thetaS3, phiS3)
return BrS3IMF, BtS3IMF, BpS3IMF
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure CheckIfInsideMappedMp
#### ;--------------------------------------------------------------------------------------
#### pro CheckIfInsideMappedMp, ctime, XS3, YS3, ZS3, ZNS3, ianswer
[docs]def CheckIfInsideMappedMp(ctime, XS3, YS3, ZS3, ZNS3):
#### vecin = (vecou = dblarr(3))
vecin = np.zeros([3])
vecou = np.zeros([3])
####
#### vecin = [XS3, YS3, ZS3 + ZNS3]
vecin = [XS3, YS3, ZS3 + ZNS3]
#### JROT, 's3c', 'jsm', vecin, vecou, ctime
vecou = jrot('s3c', 'jsm', vecin, ctime)
#### XJSOold = vecou[0]
#### YJSOold = vecou[1]
#### ZJSOold = vecou[2]
XJSOold = vecou[0]
YJSOold = vecou[1]
ZJSOold = vecou[2]
#### ;print, ZJSOold
#### ;stop
ianswer = checkIfInsideMagnetop(XJSOold, YJSOold, ZJSOold)
return ianswer
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure checkIfInsideMagnetop
#### ;
#### ; created by Mariel Desroche 11-17-2008
#### ; checks whether the position x0, y0, z0 is in the magnetopause
#### ; defined by Joy et al 2002 polynomial
#### ; Pd is the solar wind dynamic pressure (nPa)
#### ; returns 0 = outside mp, 1 = inside mp
#### ;--------------------------------------------------------------------------------------
#### ;pro checkIfInsideMagnetop, x, y, z, ianswer
#### ; ; initialize answer
#### ; ianswer = 0
####
#### ; Pd = 0.04d
#### ; Pd4 = 1d/sqrt(sqrt(Pd))
#### ; ; Magnetopause Parameters
#### ; A = -0.134d + 0.488d*Pd4
#### ; B = -0.581d - 0.225d*Pd4
#### ; C = -0.186d - 0.016d*Pd4
#### ; D = -0.014d + 0.096d*Pd
#### ; E = -0.814d - 0.811d*Pd
#### ; F = -0.050d + 0.186d*Pd
####
#### ; x0 = x/120d
#### ; y0 = y/120d
#### ; z0 = z/120d
####
#### ; aa = -0.186d - 0.016d*Pd4
#### ; bb = -0.581d - 0.225d*Pd4
#### ; cc = -0.134d + 0.488d*Pd4 - z0*z0
####
#### ; x_max = (-bb - sqrt(bb*bb - 4*aa*cc))/(2*aa)
####
#### ; aa = -0.814d - 0.811d*Pd
#### ; bb = -0.014d + 0.096d*Pd + (-0.050d + 0.186d*Pd)*x0
#### ; cc = (-0.186d - 0.016d*Pd4)*x0*x0 + (B = -0.581d - 0.225d*Pd4)*x0 - 0.134d + 0.488d*Pd4 - z0*z0
####
#### ; quad = bb*bb - 4*aa*cc
#### ; y_pos = (-bb + sqrt(quad))/(2*aa)
#### ; y_neg = (-bb - sqrt(quad))/(2*aa)
####
#### ; if (x0 lt x_max) && (y0 gt y_pos) && (y0 lt y_neg) then ianswer = 1
#### ;end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure checkIfInsideMagnetop
#### ;
#### ; original magnetopause check program... still being tested and under construction
#### ;--------------------------------------------------------------------------------------
#### pro checkIfInsideMagnetop, x, y, z, answer
[docs]def checkIfInsideMagnetop(x, y, z):
#### ; init answer
#### answer = 1
answer = 1
####
#### ; Symmetric Magnetopause Parameter
#### A = 13779.1546267827d
#### B = -130.0469911647737d
#### C = -0.2216972473764874d
#### D = 0d
#### E = -0.8453162863101464d
#### F = 0d
A = 13779.1546267827
B = -130.0469911647737
C = -0.2216972473764874
D = 0.
E = -0.8453162863101464
F = 0.
####
#### phi = atan(z,y)
#### rho_in = sqrt(z*z + y*y)
#### cos_phi = cos(phi)
#### sin_phi = sin(phi)
phi = np.arctan2(z,y)
rho_in = np.sqrt(z*z + y*y)
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
####
#### ; check if we are on the dayside
#### if (x ge 0d) then begin
if x >= 0.:
#### aa = E*cos_phi*cos_phi - sin_phi*sin_phi + F*cos_phi*sin_phi
aa = E*cos_phi*cos_phi - sin_phi*sin_phi + F*cos_phi*sin_phi
#### bb = D*cos_phi
bb = D*cos_phi
#### cc = A + B*x + c*x*x
cc = A + B*x + C*x*x
#### term = bb*bb - 4*aa*cc
term = bb*bb - 4*aa*cc
#### sqrt_term = sqrt(term)
sqrt_term = np.sqrt(term)
#### rho1 = (-bb + sqrt_term)/(2d*aa)
rho1 = (-bb + sqrt_term)/(2.*aa)
#### rho2 = (-bb - sqrt_term)/(2d*aa)
rho2 = (-bb - sqrt_term)/(2.*aa)
#### ; The right solution has rho positive (i.e. the bigger of the two)
#### rho = rho1
rho = rho1
#### if (rho2 gt rho1) then rho = rho2
if (rho2 > rho1): rho = rho2
#### if (term lt 0d) || (rho lt 0d) || (rho_in gt rho) then answer = 0
if (term < 0.) or (rho < 0.) or (rho_in > rho): answer = 0
#### endif else begin
else:
#### ; Now the nightside
#### aa = E*cos_phi*cos_phi - sin_phi*sin_phi + F*cos_phi*sin_phi
aa = E*cos_phi*cos_phi - sin_phi*sin_phi + F*cos_phi*sin_phi
#### bb = D*cos_phi
bb = D*cos_phi
#### cc = A + B*x + c*x*x
cc = A + B*x + C*x*x
#### term = bb*bb - 4*aa*cc
term = bb*bb - 4*aa*cc
#### sqrt_term = sqrt(term)
sqrt_term = sqrt(term)
#### rho1 = (-bb + sqrt_term)/(2d*aa)
rho1 = (-bb + sqrt_term)/(2.*aa)
#### rho2 = (-bb - sqrt_term)/(2d*aa)
rho2 = (-bb - sqrt_term)/(2.*aa)
#### ; The right solution has rho positive (i.e. the bigger of the two)
#### rho = rho1
rho = rho1
#### if (rho2 gt rho1) then rho = rho2
if (rho2 > rho1): rho = rho2
#### if (term lt 0d) || (rho lt 0d) || (rho_in gt rho) then answer = 0
if (term < 0.) or (rho < 0.) or (rho_in > rho): answer = 0
#### endelse
return answer
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure tail_mag_notilt_allmodes
#### ;
#### ; Modified from csheet_deform to calculate the field of only one point
#### ; The subroutine first calculates the Brho and Bz components at the stretched location
#### ; and then multiplies with the stretch matrices to calculate Bphi
#### ; alat,colat RLT are in radians, RLT is measured from Noon
#### ;--------------------------------------------------------------------------------------
#### pro tail_mag_notilt_all_modes, rho, phi, zz, RLT, Brcs, Bpcs, Bzcs, Mode
[docs]def tail_mag_notilt_all_modes(rho, phi, zz, RLT, Mode):
'''
According to the IDL code, the unit test is generated.
>>> rho = 29.774504
>>> phi = -2.6928876
>>> zz = 5.7614139
>>> RLT = 2.7842312
>>> Mode = 7
>>> Brcs, Bpcs, Bzcs = tail_mag_notilt_all_modes(rho, phi, zz, RLT, Mode)
>>> print('%.6f' % Brcs)
20.816283
>>> print('%.6f' % Bpcs)
-3.551422
>>> print('%.6f' % Bzcs)
15.714019
'''
_logger.debug('tail_mag_notilt_all_modes: rho, phi, zz, RLT, Mode =')
_logger.debug('{} {} {} {} {}'.format(rho, phi, zz, RLT, Mode))
#### f = (beta = dblarr(6,6))
f = np.zeros([6, 6])
beta = np.zeros([6, 6])
#### X = dblarr(6)
X = np.zeros([6])
#### isumno = 6
isumno = 6
####
#### ; begin Hannes Modes
#### C = [-28.402612d, -18.160631d, 14.543719d, 12.873892d, -5.590983d, 3.461748d] ; short model
C = [-28.402612, -18.160631, 14.543719, 12.873892, -5.590983, 3.461748] # short model
####
#### ; ring current
#### f[0,*] = [188.880535100774d, 170.216250249208d, -232.209387453753d, 17729.2564306361d, -17426.2836367108d, -176.740242190339d]
f[0,:] = [188.880535100774, 170.216250249208, -232.209387453753, 17729.2564306361, -17426.2836367108, -176.740242190339]
#### beta[0,*] = [59.2554626464844d, 14.3016757965088d, 3.41575503349304d, 5.47725486755371d, 5.55856323242188d, 30.7563076019287d]
beta[0,:] = [59.2554626464844, 14.3016757965088, 3.41575503349304, 5.47725486755371, 5.55856323242188, 30.7563076019287]
#### ; L = 6
#### f[1,*] = [34.5422110518d, 6.0536537495d, 414.2147150061d, -155.3647545840d, 316.3317439470d, -73.2913775361d]
f[1,:] = [34.5422110518, 6.0536537495, 414.2147150061, -155.3647545840, 316.3317439470, -73.2913775361]
#### beta[1,*] = [59.2781028748d, 16.4449157715d, 4.1893663406d, 3.0382030010d, 9.1110754013d, 30.9898948669d]
beta[1,:] = [59.2781028748, 16.4449157715, 4.1893663406, 3.0382030010, 9.1110754013, 30.9898948669]
#### ; L = 12
#### f[2,*] = [-165.1020233300d, 554.5217109532d, 315.5180924885d, -125.0952135958d, 449.1102649649d, 371.9172713735d]
f[2,:] = [-165.1020233300, 554.5217109532, 315.5180924885, -125.0952135958, 449.1102649649, 371.9172713735]
#### beta[2,*] = [55.9542961121d, 14.5898332596d, 3.9975905418d, 2.8588902950d, 8.6061649323d, 21.9151973724d]
beta[2,:] = [55.9542961121, 14.5898332596, 3.9975905418, 2.8588902950, 8.6061649323, 21.9151973724]
#### ; L = 24
#### f[3,*] = [591.9706482108d, 1376.0166168630d, 298.8401265534d, -156.9445820815d, 521.9904943354d, 2387.3415724416d]
f[3,:] = [591.9706482108, 1376.0166168630, 298.8401265534, -156.9445820815, 521.9904943354, 2387.3415724416]
#### beta[3,*] = [59.0788803101d, 16.2336521149d, 3.5862164497d, 2.7877552509d, 8.4433155060d, 30.6951904297d]
beta[3,:] = [59.0788803101, 16.2336521149, 3.5862164497, 2.7877552509, 8.4433155060, 30.6951904297]
#### ; L = 48
#### f[4,*] = [9180.3492664995d, 1679.4592666407d, 274.9968669893d, -144.4824637207d, 532.5683191252d, 4540.5253277291d]
f[4,:] = [9180.3492664995, 1679.4592666407, 274.9968669893, -144.4824637207, 532.5683191252, 4540.5253277291]
#### beta[4,*] = [59.2788314819d, 16.4751930237d, 3.5809681416d, 2.7559804916d, 8.5185470581d, 31.0038948059d]
beta[4,:] = [59.2788314819, 16.4751930237, 3.5809681416, 2.7559804916, 8.5185470581, 31.0038948059]
#### ; L = 96
#### f[5,*] = [28927.9037953983d, 2214.0596352902d, 242.5405186300d, -108.6915959578d, 609.8932996003d, 7437.5789048160d]
f[5,:] = [28927.9037953983, 2214.0596352902, 242.5405186300, -108.6915959578, 609.8932996003, 7437.5789048160]
#### beta[5,*] = [72.2351608276d, 18.0498085022d, 3.7080323696d, 2.6697711945d, 8.9703912735d, 35.2125587463d]
beta[5,:] = [72.2351608276, 18.0498085022, 3.7080323696, 2.6697711945, 8.9703912735, 35.2125587463]
#### ; end Hannes Modes
####
#### ;! P = 3.0d-3
#### ;! P = 0d
#### P = 5.73d-3
P = 5.73e-3
####
#### if (Mode gt 6) || (Mode lt 1) then begin
#### istartLoop = 0
#### istopLoop = 5
if Mode > 6 or Mode < 1:
istartLoop = 0
istopLoop = 5
#### endif else begin
#### istartLoop = Mode - 1
#### istopLoop = istartLoop
else:
istartLoop = Mode - 1
istopLoop = istartLoop
#### endelse
####
#### ; We work in the dipole magnetic coordinate system
#### ; In this subroutine:
#### ; index I = 1,7 corresponds to each component of a mode
#### ; index L = 1,6 denotes the six L modes
#### ; index j sums over observations (Br and Btheta) and is twice the number of observations
#### ; index M is used to obtain least squares equations (same dimension as L)
#### PI = 3.14159265358979d
PI = 3.14159265358979
#### PI2 = PI/2d
PI2 = PI/2.
#### DEGREE = 180d/PI
DEGREE = 180./PI
#### RADIAN = PI/180d
RADIAN = PI/180.
#### drho = 0.05d
drho = 0.05
#### dz = 0.05d
dz = 0.05
#### phi = phi
phi = phi
#### ; D is half thickness of current sheet
#### ;! D = 2.0
#### D0 = 11.0d
D0 = 11.0
#### D1 = 9.0d
D1 = 9.0
#### D = D0 + D1*cos(RLT - 1.5d*PI)
D = D0 + D1*np.cos(RLT - 1.5*PI)
#### D = 2.0d
D = 2.0
####
#### RMSBr = 0d
#### RMSBt = 0d
#### RMSBp = 0d
RMSBr = 0.
RMSBt = 0.
RMSBp = 0.
#### ; Start calculating the field
#### Brhodm = 0d
#### Bzdm = 0d
#### Bpdm = 0d
Brhodm = 0.
Bzdm = 0.
Bpdm = 0.
####
#### rhomag = rho
#### ZMAG = zz
#### Z = ZMAG
rhomag = rho
ZMAG = zz
Z = ZMAG
####
#### ; First calculate Brhod
#### ;! Do 2 L = mode,mode
#### for L = istartLoop, istopLoop do begin
for L in range(istartLoop, istopLoop + 1): # FOR python is exclusive end, but inclusive for IDL.
#### ZM = abs(Z - dz)
ZM = abs(Z - dz)
#### if (ZM lt D) then ZM = 0.5*(ZM*ZM/D + D)
if ZM < D:
ZM = 0.5 * (ZM * ZM / D + D)
#### ZP = abs(Z + dz)
ZP = abs(Z + dz)
#### if (ZP lt D) then ZP = 0.5*(ZP*ZP/D + D)
if ZP < D:
ZP = 0.5 * (ZP * ZP / D + D)
#### xlpp = 0d
#### xlpm = 0d
xlpp = 0.
xlpm = 0.
#### for i = 0, isumno - 1 do begin
for i in range(0, isumno):
#### S1p = sqrt( (beta(L,i) + ZP)*(beta(L,i) + ZP) + (RHOMAG + Beta(L,i))*(RHOMAG + Beta(L,i)) )
#### S2p = sqrt( (beta(L,i) + ZP)*(beta(L,i) + ZP) + (RHOMAG - Beta(L,i))*(RHOMAG - Beta(L,i)) )
#### S1m = sqrt( (beta(L,i) + ZM)*(beta(L,i) + ZM) + (RHOMAG + Beta(L,i))*(RHOMAG + Beta(L,i)) )
#### S2m = sqrt( (beta(L,i) + ZM)*(beta(L,i) + ZM) + (RHOMAG - Beta(L,i))*(RHOMAG - Beta(L,i)) )
S1p = np.sqrt( (beta[L,i] + ZP)*(beta[L,i] + ZP) + (rhomag + beta[L,i])*(rhomag + beta[L,i]) )
S2p = np.sqrt( (beta[L,i] + ZP)*(beta[L,i] + ZP) + (rhomag - beta[L,i])*(rhomag - beta[L,i]) )
S1m = np.sqrt( (beta[L,i] + ZM)*(beta[L,i] + ZM) + (rhomag + beta[L,i])*(rhomag + beta[L,i]) )
S2m = np.sqrt( (beta[L,i] + ZM)*(beta[L,i] + ZM) + (rhomag - beta[L,i])*(rhomag - beta[L,i]) )
#### tp = 2d*Beta(L,i)/(S1p + S2p)
#### tm = 2d*Beta(L,i)/(S1m + S2m)
tp = 2.*beta[L,i]/(S1p + S2p)
tm = 2.*beta[L,i]/(S1m + S2m)
#### AAp = tp*sqrt(1d - tp*tp)/(S1p*S2p)
#### AAm = tm*sqrt(1d - tm*tm)/(S1m*S2m)
AAp = tp*np.sqrt(1. - tp*tp)/(S1p*S2p)
AAm = tm*np.sqrt(1. - tm*tm)/(S1m*S2m)
#### xlpp = xlpp + f(L,i)*AAp*rhomag
#### xlpm = xlpm + f(L,i)*AAm*rhomag
xlpp = xlpp + f[L,i]*AAp*rhomag
xlpm = xlpm + f[L,i]*AAm*rhomag
#### endfor
#### dxpldz = (xlpp - xlpm)/(2d*dz)
dxpldz = (xlpp - xlpm)/(2. * dz)
#### X[L] = -dxpldz
X[L] = -dxpldz
#### Brhodm = Brhodm + X[L]*C[L]
Brhodm = Brhodm + X[L]*C[L]
#### endfor
####
#### ; Now calculate BZD
#### ; Now corrected for the current sheet derivative term. 11/2002.
#### ;! Do 20 L = Mode,Mode
#### for L = istartLoop, istopLoop do begin
for L in range(istartLoop, istopLoop + 1):
#### rhom = RHOMAG - drho
#### rhop = RHOMAG + drho
rhom = rhomag - drho
rhop = rhomag + drho
####
#### ; First calculate it for rhom
#### rhos3m = rho - drho
rhos3m = rho - drho
#### ZNM = 0d
ZNM = 0.
#### ZM = abs(ZMAG - ZNM)
ZM = abs(ZMAG - ZNM)
#### if (ZM lt D) then ZM = 0.5d*(ZM*ZM/D + D)
if ZM < D: ZM = 0.5*(ZM*ZM/D + D)
####
#### ; Now calculate it for rhop
#### rhos3p = rho + drho
rhos3p = rho + drho
#### ZNP = 0d
ZNP = 0.
#### ZP = abs(ZMAG - ZNP)
ZP = abs(ZMAG - ZNP)
#### if (ZP lt D) then ZP = 0.5d*(ZP*ZP/D + D)
if ZP < D : ZP = 0.5*(ZP*ZP/D + D)
#### xlpp = 0d
#### xlpm = 0d
xlpp = 0.
xlpm = 0.
#### for i = 0, isumno - 1 do begin
for i in range(0, isumno):
#### S1p = sqrt( (beta(L,i) + ZP)*(beta(L,i) + ZP) + (rhop + Beta(L,i))*(rhop + Beta(L,i)) )
#### S2p = sqrt( (beta(L,i) + ZP)*(beta(L,i) + ZP) + (rhop - Beta(L,i))*(rhop - Beta(L,i)) )
#### S1m = sqrt( (beta(L,i) + ZM)*(beta(L,i) + ZM) + (rhom + Beta(L,i))*(rhom + Beta(L,i)) )
#### S2m = sqrt( (beta(L,i) + ZM)*(beta(L,i) + ZM) + (rhom - Beta(L,i))*(rhom - Beta(L,i)) )
S1p = np.sqrt( (beta[L,i] + ZP)*(beta[L,i] + ZP) + (rhop + beta[L,i])*(rhop + beta[L,i]) )
S2p = np.sqrt( (beta[L,i] + ZP)*(beta[L,i] + ZP) + (rhop - beta[L,i])*(rhop - beta[L,i]) )
S1m = np.sqrt( (beta[L,i] + ZM)*(beta[L,i] + ZM) + (rhom + beta[L,i])*(rhom + beta[L,i]) )
S2m = np.sqrt( (beta[L,i] + ZM)*(beta[L,i] + ZM) + (rhom - beta[L,i])*(rhom - beta[L,i]) )
#### tp = 2d*Beta(L,i)/(S1p + S2p)
#### tm = 2d*Beta(L,i)/(S1m + S2m)
tp = 2.*beta[L,i]/(S1p + S2p)
tm = 2.*beta[L,i]/(S1m + S2m)
#### AAp = tp*sqrt(1d - tp*tp)/(S1p*S2p)
#### AAm = tm*sqrt(1d - tm*tm)/(S1m*S2m)
AAp = tp*np.sqrt(1. - tp*tp)/(S1p*S2p)
AAm = tm*np.sqrt(1. - tm*tm)/(S1m*S2m)
#### xlpp = xlpp + f(L,i)*AAp*rhop
#### xlpm = xlpm + f(L,i)*AAm*rhom
xlpp = xlpp + f[L,i]*AAp*rhop
xlpm = xlpm + f[L,i]*AAm*rhom
#### endfor
#### dxpldr = (rhop*xlpp - rhom*xlpm)/(2d*drho)
dxpldr = (rhop*xlpp - rhom*xlpm)/(2.*drho)
#### X[L] = dxpldr/RHOMAG
X[L] = dxpldr/rhomag
#### Bzdm = Bzdm + X[L]*C[L]
Bzdm = Bzdm + X[L]*C[L]
#### endfor
#### Bpdm = 0d
Bpdm = 0.
#### ; Now add the stretch field to Bpdm
#### Bpdm = Bpdm - p*rhomag*Brhodm
Bpdm = Bpdm - P*rhomag*Brhodm
#### ;! PH_MAG=PH_MAG-p*rhomag
####
#### ; Now rotate into system III coordinates.
#### ;! Bphi_N=0 !changed Hannes
#### Brcs = Brhodm
#### Bpcs = Bpdm
#### Bzcs = Bzdm
Brcs = Brhodm
Bpcs = Bpdm
Bzcs = Bzdm
#### ;! Bpdifm=BPHM
#### ;! Btdifm=BTHM
_logger.debug('tail_mag_notilt_all_modes: Brcs, Bpcs, Bzcs=')
_logger.debug('{} {} {}'.format(Brcs, Bpcs, Bzcs))
return Brcs, Bpcs, Bzcs
#### end
####
#### ;--------------------------------------------------------------------------------------
#### ; procedure JOVIAN_VIP4_no_dipole
#### ;
#### ; created by K.K.Khurana 1996
#### ; BASED ON THE SUBROUTINE IGRF WRITTEN BY N.A. TSYGANENKO (1979)
#### ; THE LEFT HANDED COODINATE SYSTEM SYSTEM III IS CHANGED TO A RIGHT
#### ; HANDED SYSTEM BY FEEDING IN -F.
#### ;
#### ; INPUT: NM (INTEGER)- MAXIMUM ORDER OF HARMONICS TAKEN
#### ; INTO ACCOUNT (NM.LE.12)
#### ;
#### ; R,T,F (REAL)- POSITION OF DESIRED FIELD VALUE IN
#### ; SPHERICAL GEOGRAPHIC COORDINATE SYSTEM
#### ; (R IN EARTH RADII, COLATITUDE T AND LONGITUDE F IN RADIANS)
#### ;
#### ; OUTPUT: BR,BT,BF (REAL)- COMPONENTS OF THE INTERNAL PORTION
#### ; OF THE MAIN MAGNETIC FIELD IN
#### ; SPHERICAL GEOGRAPHICAL COORD SYSTEM (VALUES GIVEN IN GAMMA)
#### ;
#### ; CALCULATES COMPONENTS OF MAIN JOVIAN FIELD IN SPHERICAL SYSTEM III (1965)
#### ; COORD SYSTEM BASING ON SPHERICAL HARMONIC COEFFICIENTS GIVEN BY ACUNA
#### ; ET.AL. (1983) IN DESSLER'S BOOK.
#### ; Now changed to O6 model. April, 1996.
#### ; MAXIMUM ORDER OF HARMONICS TAKEN INTO ACCOUNT (NOT MORE THAN 03)
#### ; R,T,F ARE SPHERICAL COORDS OF THE POINT (R IN RJ,COLATITUDE T AND
#### ; LONGITUDE F IN RADIANS), FIELD COMPONENTS BR,BT,BF
#### ; ARE CALCULATED IN NANOTESLA(= 1 GAMMA)
#### ;--------------------------------------------------------------------------------------
#### pro JOVIAN_VIP4_no_dipole, NM, r, theta, phi, BR, BT, BF
[docs]def JOVIAN_VIP4_no_dipole(NM, r, theta, phi):
#### REC = replicate(1d, 91)
REC = np.ones([91])
#### A = (B = dblarr(13))
A = np.zeros([13])
B = np.zeros([13])
#### d76 = dblarr(76)
d76 = np.zeros([76])
####
#### ; VIP4 no dipole component
#### G = [0d, 0d, 0d, -5100d, -61900d, 49700d, -1600d, -52000d, 24400d, $
#### -17600d, -16800d, 22200d, -6100d, -20200d, 6600d, d76]
G = np.array([0., 0., 0., -5100., -61900., 49700., -1600., -52000., 24400.,
-17600., -16800., 22200., -6100., -20200., 6600.] + d76.tolist())
#### H = [0d, 0d, 0d, 0d, -36100d, 5300d, 0d, -8800d, 40800d, -31600d, $
#### 0d, 7600d, 40400d, -16600d, 3900d, d76]
H = np.array([0., 0., 0., 0., -36100., 5300., 0., -8800., 40800., -31600.,
0., 7600., 40400., -16600., 3900.] + d76.tolist())
####
#### ;! VIP4
#### ;! G = [0.,420500.,-65900.,-5100.,-61900.,49700.,-1600.,-52000.,24400., $
#### ;! -17600.,-16800.,22200.,-6100.,-20200.,6600.,76*0.0]
#### ;! H = [0.,0.,25000.,0.,-36100.,5300.,0.,-8800.,40800.,-31600.,0.,7600., $
#### ;! 40400.,-16600.,3900,76*0.0]
#### ;! O6
#### ;! G = [0.,421800.,-66400.,-20300.,-73500.,51300.,-23300.,-7600.,16800.,-23100.,81*0.0]
#### ;! H = [0.,0.,26400.,0.,-46900.,8800.,0.,-58000.,48700.,-29400.,81*0.0]
####
#### G[0] = 0d
G[0] = 0.
#### H[0] = 0d
H[0] = 0.
#### KNM = 15
KNM = 15
#### for N = 1, 13 do begin
for N in range(1, 14):
#### N2 = (2d*N - 1d)*(2d*N - 3d)
N2 = (2. * N - 1.) * (2. * N - 3.)
#### for M = 1, N do begin
for M in range(1, N+1):
#### MN = N*(N - 1)/2 + M
MN = N * (N - 1) // 2 + M # I think this is //, not /
#### REC[MN-1] = (N - M)*(N + M - 2d)/N2
REC[MN-1] = (N - M)*(N + M - 2.)/N2
#### endfor
#### endfor
#### S = 1d
S = 1.
#### for N = 2, 13 do begin
for N in range(2, 14):
#### MN = N*(N - 1d)/2d
MN = int(N*(N - 1.)/2.)
#### S = S*(2d*N - 3d)/(N - 1d)
S = S*(2.*N - 3.)/(N - 1.)
#### G[MN] = G[MN]*S
#### H[MN] = H[MN]*S
G[MN] = G[MN]*S
H[MN] = H[MN]*S
#### P = S
P = S
#### for M = 2, N do begin
for M in range(2, N + 1):
#### AA = 1d
AA = 1.
#### if (M eq 2) then AA = 2d
if M == 2: AA = 2.
#### P = P*sqrt(AA*(N - M + 1d)/(N + M - 2d))
P = P*np.sqrt(AA*(N - M + 1.)/(N + M - 2.))
#### MNN = MN + M - 1
MNN = MN + M - 1
#### G[MNN] = G[MNN]*P
G[MNN] = G[MNN]*P
#### H[MNN] = H[MNN]*P
H[MNN] = H[MNN]*P
#### endfor
#### endfor
#### if (KNM ne NM) then begin
if KNM != NM:
#### KNM = NM
#### K = KNM + 1
KNM = NM
K = KNM + 1
#### endif
#### PP = 1d/r
PP = 1. / r
#### P = PP
P = PP
#### for N = 1, K do begin
for N in range(1, K+1):
#### P = P*PP
#### A[N-1] = P
#### B[N-1] = P*N
P = P*PP
A[N-1] = P
B[N-1] = P*N
#### endfor
#### P = 1d
P = 1.
#### D = 0d
D = 0.
#### BBR = 0d
BBR = 0.
#### BBT = 0d
BBT = 0.
#### BBF = 0d
BBF = 0.
#### cos_phi = cos(phi)
#### sin_phi = sin(phi)
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
#### cos_theta = cos(theta)
#### sin_theta = sin(theta)
cos_theta = np.cos(theta)
sin_theta = np.sin(theta)
#### for M = 1, K do begin
for M in range(1, K+1):
#### if (M eq 1) then begin
if M == 1:
#### X = 0d
#### Y = 1d
X = 0.
Y = 1.
#### endif else begin
else:
#### MM = M - 1
#### W = X
#### X = W*cos_phi + Y*sin_phi
#### Y = Y*cos_phi - W*sin_phi
MM = M - 1
W = X
X = W*cos_phi + Y*sin_phi
Y = Y*cos_phi - W*sin_phi
#### endelse
#### Q = P
Q = P
#### Z = D
Z = D
#### BI = 0d
BI = 0.
#### P2 = 0d
P2 = 0.
#### D2 = 0d
D2 = 0.
#### for N = M, K do begin
for N in range(M, K+1):
#### ;AN = A[N - 1]
#### MN = N*(N - 1)/2 + M
MN = N*(N - 1)//2 + M
#### ;E = G[MN - 1]
#### ;HH = H[MN - 1]
#### W = G[MN - 1]*Y + H[MN - 1]*X
W = G[MN - 1]*Y + H[MN - 1]*X
#### ;if (ABS(P2) lt 1d-38) then P2 = 0d
#### ;if (ABS(Q) lt 1d-38) then Q = 0d
#### BBR = BBR + B[N-1]*W*Q
BBR = BBR + B[N-1]*W*Q
#### BBT = BBT - A[N - 1]*W*Z
BBT = BBT - A[N - 1]*W*Z
#### if (M ne 1) then begin
if M != 1:
#### QQ = Q
#### if (sin_theta lt 1d-5) then QQ = Z
#### BI = BI + A[N - 1]*(G[MN - 1]*X - H[MN - 1]*Y)*QQ
QQ = Q
if sin_theta < 1e-5: QQ = Z
BI = BI + A[N - 1]*(G[MN - 1]*X - H[MN - 1]*Y)*QQ
#### endif
#### XK = REC[MN-1]
XK = REC[MN - 1]
#### DP = cos_theta*Z - sin_theta*Q - XK*D2
DP = cos_theta*Z - sin_theta*Q - XK*D2
#### PM = cos_theta*Q - XK*P2
PM = cos_theta*Q - XK*P2
#### D2 = Z
D2 = Z
#### P2 = Q
P2 = Q
#### Z = DP
Z = DP
#### Q = PM
Q = PM
#### endfor
#### D = sin_theta*D + cos_theta*P
D = sin_theta*D + cos_theta*P
#### P = sin_theta*P
P = sin_theta*P
#### if (M ne 1) then begin
if M != 1:
#### BI = BI*MM
#### BBF = BBF + BI
BI = BI*MM
BBF = BBF + BI
#### endif
#### endfor
#### BR = BBR
BR = BBR
#### BT = BBT
BT = BBT
#### if (sin_theta gt 1d-5) then begin
if sin_theta > 1e-5:
#### BF = BBF/sin_theta
BF = BBF/sin_theta
#### endif else begin
else:
#### if (cos_theta lt 0d) then BBF = -BBF
#### BF = BBF
if cos_theta < 0: BBF = -BBF
BF = BBF
#### endelse
return (BR, BT, BF)
#### end
####
####
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')