JuliaSpace / SatelliteToolbox.jl

A toolbox for satellite analysis written in julia language.
MIT License
249 stars 33 forks source link

Complex Values Causing Issues in JB2008 #39

Closed abukowski21 closed 4 years ago

abukowski21 commented 4 years ago

Something strange is happening with complex variables in the JB2008 model. I tried changing all the erroneous exponentials to Complex(x)^y, but that didn't fix it. I'm trying to compare GOCE neutral density data to JB2008's predictions in May 2010. This error seems to happen for most of that month (~180,000/~242,000 times). I haven't seen any errors like this for several other runs.

ERROR: DomainError with -0.9862171602726729: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar. Stacktrace: [1] throw_exp_domainerror(::Float64) at ./math.jl:36 [2] ^ at ./math.jl:849 [inlined] [3] jb2008(::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64) at /Users/abukowski/.julia/packages/SatelliteToolbox/blx7k/src/earth/atmospheric_models/jb2008/jb2008.jl:201

jb2008(DatetoJD(2010, 5, 1, 0, 44, 18), deg2rad(355.94), deg2rad(83.40), 292000.00, 80.80,76.40,78.70,73.10,79.00,75.90,86.80,79.60,31.00)

Any ideas?

ronisbr commented 4 years ago

The second parameter of jb2008 is the latitude and the third is the longitude. You can see the documentation by typing: ?jb2008:

help?> jb2008
search: jb2008 JB2008_Output

  function jb2008(JD::Number, glat::Number, glon::Number, h::Number)
  function jb2008(JD::Number, glat::Number, glon::Number, h::Number, F10::Number, F10ₐ::Number, S10::Number, S10ₐ::Number, M10::Number, M10ₐ::Number, Y10::Number, Y10ₐ::Number, DstΔTc::Number)

  Compute the atmospheric density using the Jacchia-Bowman 2008 (JB2008) model.

  If the space indices are not provided (first call), then they will be obtained from the online database. In this case,
  the function init_space_indices() must be called first and the function will throw an exception if the selected JD is
  outside of the available data.

  This model is a product of the Space Environment Technologies, more information can be seen in the websites:

  http://sol.spacenvironment.net/jb2006/

  http://sol.spacenvironment.net/jb2008/

  Args
  ≡≡≡≡≡≡

    •    JD: Julian day.

    •    glat: Geocentric latitude [rad].

    •    glon: Geocentric longitude [rad].

    •    h: Altitude [m].

    •    F10: 10.7-cm solar flux [10⁻²² W/(M² Hz)] (Tabular time 1 day earlier).

    •    F10ₐ: 10.7-cm averaged solar flux, 81-day centered on input time (Tabular time 1 day earlier).

    •    S10: EUV index (26-34 nm) scaled to F10.7 (Tabular time 1 day earlier).

    •    S10ₐ: EUV 81-day averaged centered index (Tabular time 1 day earlier).

    •    M10: MG2 index scaled to F10.7 (Tabular time 2 days earlier).

    •    M10ₐ: MG2 81-day averaged centered index (Tabular time 2 days earlier).

    •    Y10: Solar X-ray & Lya index scaled to F10.7 (Tabular time 5 days earlier).

    •    Y10ₐ: Solar X-ray & Lya 81-day averaged centered index (Tabular time 5 days earlier).

    •    DstΔTc: Temperature variation related to the Dst.

  Returns
  ≡≡≡≡≡≡≡≡≡

  An instance of the structure JB2008_Output with the computed values.

Changing the parameters we have:

julia> jb2008(DatetoJD(2010, 5, 1, 0, 44, 18), deg2rad(83.40), deg2rad(355.94), 292000.00, 80.80,76.40,78.70,73.10,79.00,75.90,86.80,79.60,31.00)
JB2008_Output{Float64}
  nN2: Float64 4.999327095536922e13
  nO2: Float64 2.581844106708866e12
  nO: Float64 3.982511703541514e14
  nAr: Float64 8.795664543367258e9
  nHe: Float64 6.52909052144658e12
  nH: Float64 1.9751806978697314e11
  rho: Float64 1.3086707143236572e-11
  T_exo: Float64 784.9497202580523
  Tz: Float64 770.6360757426783