[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 [1], GetColors [2], GetColors [3], 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:
NewcombSun
PURPOSE:
Calculates the true ecliptic coordinates (longitude, latitude and
distance) of the Sun according to Newcomb's theory of solar motion.
See O. Montenbruck, Practical Ephemeris Calculations, par. 4.1.2, p. 66.
CATEGORY:
smei/gen/idl/ephem
CALLING SEQUENCE:
FUNCTION NewcombSun, UT, $
Precision = Precision , $
Longitude = Longitude , $
Latitude = Latitude , $
Distance = Distance , $
diskradius = DiskRadius , $
parallax = Parallax , $
degrees = Degrees , $
heliocentric=Heliocentric
INPUTS:
UT array; type: standard time structure
times (UT)
OPTIONAL INPUTS:
/Precision if set, perturbations by several planets are
taken into account
/Longitude if set, return ecl. longitude
/Latitude if set, return ecl. latitude
/Distance if set, return Sun-Earth distance
(in neither of the tree keywords LONGITUDE,LATITUDE,
and DISTANCE is set, all three quantities are returned)
/Degrees if set, return angles in degrees (default: radians)
/Heliocentric if set, the longitude and latitude of Earth are returned
(180 degrees is added to the longitude, the latitude
is multiplied by -1, the distance is the same)
OUTPUTS:
Result array; type: double
cliptic coordinates of the Sun
ecliptic longitude and latitude (degrees or radians)
and Sun-Earth distance in AU
If only one coordinate (longitude, latitude or distance) is
requested then the structure of the output array the same as for T.
If more than one coordinate is requested the output array has an
additional leading dimension with 2 or 3 elements. The leading dimension
contains longitude, latitude and distance, in that order.
OPTIONAL OUTPUTS:
Diskradius radius of Sun in arcsec
Parallax horizontal parallax in arcsec
INCLUDE:
@compile_opt.pro ; On error, return to caller
CALLS: ***
InitVar, SyncDims, TimeGet, ToRadians
CALLED BY:
GetColors [2], GetColors [3], PlotEarthSkymap [2], PlotEarthSkymap [3]
PlotEarthSkymap [4], SkyDriveIn, ipsv [1], ipsv [2], ipsv [3], jpl_test, lospos
test_td, view3d, vu_thomson_antifish, vu_thomson_hammer
SEE ALSO:
jpl_eph, jpl_test
RESTRICTIONS:
The calculation of Diskradius and Parallax needs the Sun-Earth
distance. Hence the keyword /Distance must be set or all three
keywords /Distance, /Longitude and /Latitude should be omitted.
PROCEDURE:
DLP long period perturbation of solar mean longitude and mean
anomaly
G mean anomaly of the Sun
LO = DLO+SLO/3600 = mean longitude of Sun
DL difference between true and mean solar longitude according to
two body problem
The precision calculation includes perturbations by planets, oscillatory
terms with amplitudes less than 8".
The alternative way to calculate the geocentric location of the
Sun or the heliocentric location of the Earth is to use jpl_eph
(or equivalently, big_eph). In most of our IDL code NewcombSun
calls, have been replaced by big_eph calls.
The following two calls should give essentially the same result:
p = NewcombSun(tt,/precision,/degrees)
p = big_eph(tt,body='sun',center='earth',/degrees,/to_sphere,/to_ecliptic,/precess)
A check with the Astronomical Almanac shows:
(All times are at midnight on the day given)
Normal Precision Almanac
1971 Jan 1 Lng 279.9114 279.91575 279 54' 57.0" (279.91583)
Lat 0. 0.38" 0.42"
Dist 0.9832947 0.9833006 0.9832998
1974 Jan 1 Lng 280.1882 280.1874 280 11' 14.94" (280.1875)
Lat 0.
Dist
1979 Jan 1 Lng 279.9702 279.9708 279 58' 14.90" (279.9708)
Lat 0.
Dist
1989 Jan 1 Lng 280.5535 280.5525 280 33' 10.44" (280.5529)
Lat 0.
Dist
MODIFICATION HISTORY:
1989, Paul Hick (MPAE,UCSD)
10/24/91, Tom Davidson : accuracy test
AUG-1993, Paul Hick, extension of SunEclLong
FEB-1998, Paul Hick (UCSD/CASS, pphick@ucsd.edu)
modified to accept multi-dimensional input arrays
added keyword /heliocentric
replaced keyword /radians by keyword /degrees
[Previous]
[Next]
NAME:
next_contiguous_group
CALLING SEQUENCE:
FUNCTION next_contiguous_group, indx, count=count
INCLUDE:
@compile_opt.pro
COMMON:
common next_contiguous_group_common
n
m
p
nmin
CALLS: ***
InitVar, destroyvar, init_contiguous_group
CALLED BY:
init_contiguous_group, smei_hdr_plot
[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)
[Previous]
[Next]
NAME:
nso_fe_plot
PURPOSE:
Plot synoptic maps for Fe XIV and Fe X data
CALLING SEQUENCE:
nso_fe_plot
OPTIONAL INPUTS:
All input keywords for nso_fe_read are permitted.
OPTIONAL OUTPUTS:
CALLS: ***
Carrington, InitVar, PlotSynopticMap, TimeOp, TimeSet, nso_fe_read
PROCEDURE:
> Green and red maps are calibrated using the constant calibration
for post-1989 observations.
> Temperatures are calculated by taking ratios of the green and red
arrays, i.e. this assumes that the first scan of green and red map
are taken at the same time.
MODIFICATION HISTORY:
SEP-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
[Previous]
[Next]
NAME:
nso_fe_read
PURPOSE:
Read red (Fe X), green (Fe XIV) or yellow (Ca XV) synoptic map from file
CATEGORY:
NSO
CALLING SEQUENCE:
map = nso_fe_read( /green_line, /calibrate, source=source)
INPUTS:
OPTIONAL INPUT PARAMETERS:
source=source scalar; type: string; default: $SSW_SMEI_DAT/map/sacpeak
location of Sac Peak maps
The type of data read is set using one of the following keywords:
/temperature get temperature map
/green_line get green line map (same as line = 'g')
/red_line get red line map (same as line = 'r')
/yellow_line get yellow line map (same as line = 'y')
line=line scalar; type: string
'g','r' or 'y'.
If none of these are set then line='g' is assumed.
Keyword /temperature results in two recursive calls to nso_fe_read to read a
green and a red map, and calculate the temperature map from the ratio.
The 'line' setting is also used to limit the files presented in the IDL
pickfile dialog.
rotation=rotation
scalar; type: integer; default: none
limb scalar; type: string
limb indicator, 'ww', 'we', 'ew' or 'ee'
/bestlimb if set this override the input 'limb' setting and forces
reading of the 'best limb' map.
height=height integer; type: scalar; default: 115 ('g' and 'r') or 113 ('y')
scan height in solar radii, times 100 (i.e. 115, 125, 135 or 145)
If 'rotation' and 'limb' (or /bestlimb) are set than the IDL pickfile dialog
is NOT called and a file name is constructed (also using 'line' and 'height').
/calibrate if set the green and red data are scaled using the post-1989 calibration:
Map = 3.49*0.693*Map Green
Map = 2.88*0.543*Map Red
The yellow line is never calibrated
silent=silent controls messages to screen
OUTPUTS:
map array[32,61]; type: float
line intensities for 32 daily scans (coverng a little more than
one Carrington rotation).
error scalar; type: string
error string returned by flt_read
error = '' after successfull read
line scalar; type: string
'g', 'r' or 'y'
rotation scalar; type: integer
Carrington rotation number
limb scalar; type: string
limb indicator, 'ww', 'we', 'ew' or 'ee'
height scalar; type: integer
scan height in solar radii, times 100 (i.e. 115, 125, 135 or 145)
utstart array[1]; type: time structure
start time for map (=time of first scan)
CALLS: ***
FILEPATH, GetFileSpec, InitVar, IsType, SetFileSpec, flt_read, flt_string
nso_fe_start, nso_fe_temperature
CALLED BY:
nso_fe_plot
PROCEDURE:
The default location for the Sac Peak files is directory $SSW_SMEI_DAT/map/sacpeak
MODIFICATION HISTORY:
JUL-2001, Paul Hick (UCSD/CASS)
SEP-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
Rewrite
[Previous]
[Next]
NAME:
nso_fe_start
PURPOSE:
Retrieves the exact start time for an NSO Fe XIV, Fe X or Ca XV synoptic map
CATEGORY:
gen/idl/toolbox/nso_fe
CALLING SEQUENCE:
status = nso_fe_start(rotation, limb_select, line=line, $
utstart=utstart, limb=limb, height=height)
INPUTS:
rotation scalar; type: integer
Carrington rotation number
limb_select scalar; type: string
one of 'ww','ee','ew','we'
If not set then the best limb is used (as identified by the
plus sign in the sacpeak.* files).
OPTIONAL INPUT PARAMETERS:
line=line scalar; type: string; default: 'g'
one of 'g', 'r', 'y'
coronal line id: green (Fe XIV), red (Fe X) or yellow (Ca XV)
height=height integer; type: scalar; default: 115 ('g' and 'r') or 113 ('y')
scan height in solar radii, times 100 (i.e. 115, 125, 135 or 145)
/silent controls messages
OUTPUTS:
status 0: no start time found (utstart and limb do not exist)
1: start time found
output variables utstart and limb will only exist if status=1.
OPTIONAL OUTPUT PARAMETERS:
error=error scalar; type: string
error string if status=0; null string if status=1
utstart=utstart
array[1]; type: time structure
start time of synoptic map
limb=limb scalar; type: string
one of 'ww','we','ee','ew'
Same as limb_select if it was set, or else the 'best limb'
CALLS: ***
ArrayLocation, FILEPATH, InitVar, TimeSet, destroyvar, flt_read, who_am_i
CALLED BY:
nso_fe_read
RESTRICTIONS:
The startimes are stored in three files sacpeak.xiv, sacpeak.x and sacpeak.ca
stored in the same directory ast this routine.
PROCEDURE:
MODIFICATION HISTORY:
JUL-2001, Paul Hick (UCSD/CASS)
NOV-2003, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
Complete rewrite.
[Previous]
[Next]
NAME:
nso_fe_temperature
PURPOSE:
Temperature estimate from ratio of red and green intensity
CATEGORY:
NSO
CALLING SEQUENCE:
T = nso_fe_temperature, green, red
INPUTS:
green array; type: float
green (Fe XIV) intensity
red array; type: float
red (Fe X) intensity
OUTPUTS:
T array; type: float
temperatures (in MK??)
CALLS: ***
BadValue
CALLED BY:
nso_fe_read
PROCEDURE:
Where both intensity are <= zero no temperature is
calculated (returned as !values.f_nan)
The temperature is calculated from
T = 0.138645-0.0717*alog10(red/green)
Where only the green intensity is available the temperature
is set to the maximum temperature Where only the red intensity
is available the temperature is set to the minimum temperature
MODIFICATION HISTORY:
JUL-2001, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
