[Previous]
[Next]
NAME:
nagoya_glevel
PURPOSE:
Estimate G-level from scintillation index
CATEGORY:
sat/idl/toolbox/ips
CALLING SEQUENCE:
PRO nagoya_glevel, TT, $
plotx = plotx , $
all_data = all_data , $
input = input , $
output = output , $
rnum = rnum , $
nsg = nsg , $
whole_year = whole_year , $
use_min_glevel = use_min_glevel, $
min_ppdist = min_ppdist , $
_extra = _extra
INPUTS:
TT array[1]; type: time structure
only g-levels for source observation in TT +/- 0.5 days
are updated (unless /whole_year is set).
OPTIONAL INPUT PARAMETERS:
input=input scalar; type: string; default: $SSW_SMEI_DAT/nagoya/daily
fully-qualified name of a yearly Nagoya IPS file of type
nagoya.yyyy, or directory where these files are located.
If a directory is specified then the most recent yearly file
is used.
The default is the directory that holds the yearly files
derived from the real-time Nagoya IPS data. These are the
files created and maintained by dailyips.f by ingesting
the VLIST_UCSD_* files downloaded from Nagoya daily.
/all_data if set, g-level is calculated from all data.
if not, g-level is calculated from past data.
/use_min_glevel if scintillation indices are available from multiple stations
calculated g-level for all stations and use the lowest
g-level as the final value.
/output if set, update the input file
output=output scalar; type: string; the name of a file into which the
updated data are written (instead of overwriting the input
file if /output is specified).
/plotx if set, a GIF map is created and displayed
rnum=rnum scalar; type: integer; default: 8
minimum # data point required for reliable fitting
nsg=nsg scalar; type: float; default: 1.5
criterion to eliminates highly-scattered data
from fitted line. (nsg*sigma)
/whole_year if set, all g-levels in the yearly file are recalculated
(argument TT, if set, is ignored).
min_ppdist=min_ppdist
scalar; type: float; default: 0.2
minimum point-P distance in AU. IPS observations with distances
smaller than this are not used in the calculation of the scint index
vs. point-P distance dependence of each IPS source.
OUTPUTS:
If /plotx is set then the sources in TT +/- 0.5 days are displayed in a fisheye plot.
If keyword 'output' is used then the modified G-levels are written to file.
If neither keyword is used then a 'test run' is made with only text messages
to the screen.
INCLUDE:
@compile_opt.pro ; On error, return to caller
CALLS: ***
CheckDir, FILEPATH, GetColors, InitVar, IsTime, IsType, STDDEV, STRETCH, TimeFixYear
TimeGet, TimeSet, TimeUnit, UNIQ, boost, flt_string, nagoya_plotg2d, nagoya_r_fitting
twin, txt_read
PROCEDURE:
The IPS data should be stored in yearly files with
name 'nagoya.yyyy' where yyyy is the year.
IPS DATA FORMAT: (OLD)
0 1 2 3 4 5 6 7
0123456789012345678901234567890123456789012345678901234567890123456789012345
SOURCE YRMNDY UT DIST HLA HLO GLA GLO CARR V ER G-LEVEL SC-INDX G-LEVEL SC-INDX G-LEVEL SC-INDX
1117+14 000711 7.10 0.85 6 -33 6 316 1965 303-999 0.00000 0.518E+02 0.00000 0.518E+02 0.00000 0.518E+02
The second pair g-level/scint index pair contains the g-level and scint index for Kiso.
The third pair contains the g-level and scint index for Fuji.
The first g-level/scint index pair is a 'best value' pair. Currently the Fortran program
dailyips puts the scint index for Kiso there, while the g-level gets filled
in here depending on keyword selection. By default the Kiso g-level is used;
if /use_min_glevel is set then the smallest g-level (Kiso or Fuji) is used.
Scintilation index values, S, are fitted to a power law function of the point-P distance, r:
S = 10^a * r^b
(a,b) are determined by a linear least-square fit.
The g-level, g, for an individual scintillation index follows by normalizing to the
fitted function:
g = sqrt( S/(10^a+r^b))
Note that the g-level for a fixed scintillation will change if the fit changes (i.e. if
a different set of observations is used to determine the fit.
MODIFICATION HISTORY:
JUL-2000, K. Fujiki (STELAB)
JUL-2000, Paul Hick (UCSD/CASS)
The calculation of g-levels was modified in two ways:
1. The calculation of standard deviations is now done with the IDL
routine 'stddev', rather than Fujiki's function 'sigma' (which gave
a slightly different value)
2. The sources for which the g-levels are to be recalculated are all sources
within 0.5 days of the specified time. Fujiki's original code would sometimes
miss a few sources when the time would be close to the start or end of a month.
The 'fish-eye' plot now shows the elongation in the radial direction (consistent
with the display in vu_earthskymap (Fujiki's original code used the point-P
distance as the radial coordinate).
SEP-2001, Paul Hick (UCSD/CASS)
Added the /whole_year keyword to process all observations in a yearly file in one pass.
SEP-2002, Paul Hick (UCSD/CASS)
Modified to handle files with multiple scintilation indices.
Two modifications were required: the txt_read call now uses
/test to determine the record length instead of a hardcoded length
of 76. The output section truncated records after the first
scintillation index. Now it doesn't. The current modifications
should allow additions of other scint indices in the future.
NOV-2002, Paul Hick (UCSD/CASS)
Added separate calculations of Kiso and Fuji g-levels.
Added keyword /use_min_glevel (which determines calculation of a 'best' g-level.
JUL-2003, Paul Hick (UCSD/CASS)
Fixed bug in fitting procedure for g-levels (in nagoya_r_fitting). It would crash
if only two data values are available.
APR-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
Fixed bug in TimeSet with /botime set (keywords changed some time ago)
AUG-2009, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
Added some documentation, Put two helper functions into separate files.
[Previous]
[Next]
NAME:
nagoya_plotg2d
CALLING SEQUENCE:
PRO nagoya_plotg2d, r, glevel, color=color, source=source
INCLUDE:
@compile_opt.pro ; On error, return to caller
CALLS: ***
CvPointOnLos, FishEye
CALLED BY:
nagoya_glevel
[Previous]
[Next]
NAME:
nagoya_r_fitting
CALLING SEQUENCE:
FUNCTION nagoya_r_fitting, ppdist, scindx, min_ppdist
INPUTS:
ppdist array[n]; type: float
"point-P" distance of lines of sight
scindx array[n]; type: float
scintillation index
min_ppdist scalar; type: float
only data point with point-P distance
larger than min_ppdist are used in the fit
INCLUDE:
@compile_opt.pro ; On error, return to caller
CALLS: ***
SVDFIT
CALLED BY:
nagoya_glevel
PROCEDURE:
Fits scintillation index to a power low function of the point-P
distance:
scindx = 10^a[0] * ppdist^a[1]
= 10^( a[0] + a[1]*alog10(ppdist) )
alog10(scindx) = a[0] + a[1]*alog10(ppdist)
MODIFICATION HISTORY:
AUG-2009, Paul Hick (UCSD/CASS)
Added documentation
[Previous]
[Next]
NAME:
nrSimpson
PURPOSE:
Numerical integration of function of one variable using
Simpson rule
CALLING SEQUENCE:
R = nrSimpson(F,A,B,eps=DEL,maxn=MAXN,n=N,arg=arg)
INPUTS:
F real function; declared external in caller
A,B real integration range
DEL real accuracy
MAXN integer maximum number of iterations
OUTPUTS:
R real result of integration
N integer 0: no convergence
1: if A=B (then also R=0)
CALLED BY:
IPS_WeightFnc
PROCEDURE:
> If DEL <=0 then DEL=1x10^-4 is used
> If MAXN< 2 then MAXN=2 is used
MODIFICATION HISTORY:
[Previous]
[Next]
NAME:
nrZbrac
PURPOSE:
Outward search to attempt bracketing a root of Func
CATEGORY:
Math
CALLING SEQUENCE:
OK = nrZbrac(Func, X1,X2, fixX=FixX, arg=Arg)
INPUTS:
Func string
name of function
X1,X2 array
arrays of same dimension, containing
endpoints of initial guesses for
intervals to be extended
OPTIONAL INPUTS:
fixX=FixX integer scalar 0: extend on both sides (default)
1: extend X2 side only (keeping X1 fixed)
2: extend X1 side only (keeping X2 fixed)
arg=Arg anything
passed unmodified as argument to Func
OUTPUTS:
OK byte array
same dimension as X1,X2
0: root has not been bracketed
1: root has been bracketed
X1,X2 array
updated endpoints of intervals:
Func(X1)*Func(X2) <=0 where OK=1
OPTIONAL OUTPUTS:
CALLS: ***
InitVar, SyncArgs
CALLED BY:
ThomsonMidpoint, ThomsonMidpointFar, smei_radial2theta
COMMON BLOCKS:
SIDE EFFECTS:
If no initial guess is provided i.e. X1=X2 for any interval, then
OK is set to 0 for that interval
RESTRICTIONS:
PROCEDURE:
The function Func is called as
F1 = call_function(Func, X, Arg)
with X always an array of the same dimension as X1,X2. This allows
Arg to be an array containing different input for different intervals.
The endpoints of the interval, X1 and X2, are increased by factors
of Factor=1.6 until the root is bracketed, or until nTry=50 extensions
have been performed.
MODIFICATION HISTORY:
JAN-1998, Paul Hick (UCSD/CASS)
FEB-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
Fixed bug triggered when X1 and X2 are arrays
[Previous]
[Next]
NAME:
nrZbrak
PURPOSE:
Bracket one or more zero's in a given interval for a function of one
independent variable
CATEGORY:
Math: finding roots
CALLING SEQUENCE:
nrZBrak, Func,X1,X2,N,XB1,XB2,NB, maxzero=MaxZero, arg=Arg
INPUTS:
Func string
function to be tested
X1,X2 arrays[nX]
limits of interval to be tested
N scalar; type: integer
# subintervals to be tested
OPTIONAL INPUTS:
maxzero=MaxZero
scalar; type: integer
maximum # sign changes looked for (default: 1)
arg=arg anything; passed as extra argument to Func
OUTPUTS:
NB[nX] 1-dim array
# sign changes detected
XB1[nX,MaxZero]
XB2[nX,MaxZero]
2-dim array limits of subintervals over which sign change is detected
CALLS: ***
SyncDims
PROCEDURE:
"Inward search" for roots. Subdivide a given interval [X1,X2] into
N subintervals and look for sign changes across the subintervals.
See Press et al., "Numerical Recipes", Cambridge UP (1989), par. 9.1, p. 245
[Previous]
[Next]
NAME:
nrZBrent
PURPOSE:
Returns roots of function Func between X1 and X2 (with accuracy Tol)
CATEGORY:
MATH: Numerical Recipes
CALLING SEQUENCE:
Z0 = nrZBrent(Func,X1,X2,Tol,nar=NaR,arg=Arg)
INPUTS:
Func string name of function for which roots are needed
X1,X2 arrays of same size
containing two coordinate values which
bracket the root
Tol scalar tolerance for the calculated root
OPTIONAL INPUTS:
nar=NaR scalar Not a Root: number used to indicate a
root that was not bracketed by the
input X1,X2 values
default: !Values.F_NAN
arg=Arg anything is passed unmodified to all Func calls
For instance, Arg can be an array of
the same size as X1,X2, with separate
argument values for each of the
intervals [X1,X2]
OUTPUTS:
nrZBrent array of same size as X1, X2
the calculated roots
CALLS: ***
SyncArgs
CALLED BY:
KeplerOrbit, ThomsonMidpoint, ThomsonMidpointFar, smei_radial2theta
RESTRICTIONS:
The root must be bracketed by X1 and X2, i.e. Func(X1)*Func(X2) <= 0.
The function for which the root is calculated is called in the form
F = call_function(Func, X, Arg), where X is an array of the same
size as X1 and X2.
PROCEDURE:
See Numerical Recipes, Press and Teukolsky
The main difference with the Numerical Recipes version is that
exact roots are intercepted right at the start of the iteration loop
(in Numerical Recipes they are intercepted at the convergence check).
MODIFICATION HISTORY:
JAN-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)