#82 get_mag_onera.pro comments
Some things I've found that could be tweaked/improved...
The description at the top of the file needs to be updated. Inc.
> ;* EXTMOD external field model index
> ;* 5: Olson & Pfitzer quiet [1977] (default)
>
> yet line 216 is:
> if keyword_set(EXTMOD) then kext = long(EXTMOD) else kext = 0l
> So it's setting the default to 0, not 5.
> Line 223:
> IF keyword_set(PITCH) THEN pitch = pitch ELSE pitch = [90.0d]
> I suspect this would be safer as:
> IF keyword_set(PITCH) THEN pitch = double(pitch) ELSE pitch = [90.0d]
>
>
> The procedure starts with:
> PRO get_mag_onera, tai, pos, kcompstr, kinstr, $ ;inputs
> b1, b2, b3, Bmin, Lm, Lstar, XJ, MLT, mlat, equ_fl, error, $ ;outputs
> INTMOD = INTMOD, EXTMOD = EXTMOD, MAG_DYN = MAG_DYN, $
> INT_NAME = INT_NAME, EXT_NAME = EXT_NAME, $
> SW_DECK = SW_DECK, DECK = DECK, FILEDECK = FILEDECK, $
> DO_LSTAR = DO_LSTAR, LANLSTAR = LANLSTAR, LSTAR_RES = LSTAR_RES, $
> PITCH = PITCH, VERBOSE = VERBOSE, HAVE_MDATE = HAVE_MDATE, $
> DRIFTSHELL = DRIFTSHELL, DO_IK = DO_ik, IKVALUE = ikvalue, $
> ;--cy--, 102504, make calculation of I,K(values in IKVALUE)
> TRACEMP = tracemp, MPPOS = mppos, $
> ;--cy--, 102504, trace mirror points(values in mppos)
> YEAROFFSET = YEAROFFSET, NODATA = NODATA
>
> So in order:
> INTMOD, the values is used to set kint on line 214:
> > if keyword_set(INTMOD) then kint = long(INTMOD) else kint = 0l
>
> However kint isn't used anywhere else in the code. So pointless. :-)
>
> These keywords are not used anywhere at all:
> INT_NAME, EXT_NAME, VERBOSE, HAVE_MDATE, DRIFTSHELL
>
> If DECK is set, then the code is hardwired to look in directories
> /n/projects/cda/lstar_deck
> so not useful to those outside of LANLs linux cluster.
>
> Line 222:
> > IF keyword_set(LSTAR_RES) THEN lstar_res = long(LSTAR_RES) else
> lstar_res = 0
> Safer as
> > IF keyword_set(LSTAR_RES) THEN lstar_res = long(LSTAR_RES) else
> lstar_res = 0l
> (although I think this is taken care of later, but for consistency)
>
> Similarly for NODATA, should line 214 that is currently:
> > if keyword_set(NODATA) then nodata = NODATA else nodata = -99
> be:
> > if keyword_set(NODATA) then nodata = double(NODATA) else nodata = -99.0d
>
> I'm not so sure on that last one.
Line 230 for m_in, alpha and Npa are limited to array size of 25 - which seems to be due to the onera libraries and what's hardwired in them.
Maybe 25 should just be replaced by n_elements(pitch), and/or if that's bigger than 25 then a warning message should be given to the screen.
Also when I run this I get the messages returned:
% Program caused arithmetic error: Floating divide by 0
% Program caused arithmetic error: Floating illegal operand
I've established that these errors are reported in IDL at the end of the execution of the code, and not during the code - making locating the lines that cause this hard.
I've gone through every divide sign in the IDL procedure to check it's not diving by zero - they are not.
So is it possible these errors are coming from the fortran codes called by this procedure?
Variable LStar is specified on line 364 as: Lstar = dblarr(n,25) (where n was set to 100001, another hard set for the fortran code).
If the DO_LSTAR keyword in the procedure is set to =0 this returns a big array for Lstar of size 10000x25, which seems a waste of memory when it's empty.
Line 604, pa = pitch[p] is in a for loop, for p = 0, 8 , why 8? Shouldn't that be 0,n_elements(pitch)-1 ??
And an onera bug - I suspect... if you feed it a scalar pitch angle and the field line doesn't close then you get a fill value for MLAT and LM. Fair enough. But if pitch_angle = [5, 10, 15, ... ,80, 85, 90] and the 5 degrees doesn't close then MLAT and LM are fill values. However, if pitch_angle = [90, 85, 80, ..., 15, 10, 5] (i.e. just flipped the order) then MLAT and LM are populated correctly. I assume the onera library populates LM and MLAT based on pitch_angle(0), and if it's fill it doesn't try the next pitch_angle index.
Discussion
