In this version of inscal, the native raw data format is the midas cmst structure used at Millstone Hill, and the output format is CEDAR format. The input routine, readifamrecord, can handle both the current cmst format and historical cmst and ifam file formats. In addition, the ncar graphics package is used to produce a summary plot of each analyzed record. These are stored in subdirectories of the Madrigal experiment directory and are accessible through the Madrigal Web interface.
inscal - Inscal main program
setsignals - Register exinsc to catch signals
argin - Processe input arguments
ctrlin - Read control parameters
openmhrparms - Open calibration file (e.g. mhr_parms.dat)
openmidasfile - Open input file (Cmst or ifam - any version)
cedaropen - Open Madrigal file for each record in input file
readifamrecord - Read next record of input file
postint - Postintegrate when the ifam record contains subintegrations
(e.g. median of subintegrations)
profint - smooth each lag profile
(e.g. median filter, L1 spline fit etc.)
acfit - Process single record (integration period)
for each range in record
aquik,arget - Get initial guess for fit parameters
fit1 - Fit for Ti, Te
for each iteration in least-squares fit
curfit
functn - Evaluate theroretical acf
end for
fit1 - Add CO at low altitudes or PH at high altitudes
for each iteration in least-squares fit
curfit
functn - Evaluate theroretical acf
end for
end for
ofam2cedar - put results in ofam
plotpar - Plot Ti, Te, Vo and Pol
ofam2cedar - Convert inscal native output format (ofam) to external format (Cedar)
and write to Madrigal file
end for
exinsc - Clean up, print statistics and exit
inscal argument_list
Arguments:
-cmdfile - File containing a series of arguments (commands), one
line per argument. When this argument is encountered on
the command line, cmdfile is read and its contents are
treated exactly as if they had been entered on the
command line. Typically this is the first argument in
a call to inscal is followed by a small number of
additional arguments which overide the values in
cmdfile.
-oasisdir - oasisdir where analysis results will be wriitten.
e.g.:
/usr/oasis/experiments/2003/20nov03/oasis/inscal
Default is the current oasisdir
-tag - Character string tag added to end of inscal output files.
Use #startTime to specify that tag should be the current
Unix time in seconds, e.g. "_1101217411"
-tag #startTime
-outfile - Madrigal output file. Use #tag to specify mlh_tag where
tag is specified by the tag command.
-outfile #tag
-insfile - Real-time status file. Inscal gets the input file
(otherwise specified by -infile) and output Madrigal
file (otherwise specified by -outfile) from insfile
if -runmode = "realTime".
-infile - Input file, typically a cmst file. Full path name.
-infile "infile_default'
-runmode - dataFile: Get infile and outfile from inscal arguments
realTime: Get infile and outfile from insfile
asciiIfam: Like dataFile, but infile is ascii format
-startrec - First record to analyze
-startRec 1
-endrec - Last record to analyze
-endrec 99999999
-delrec - Record increment, e.g. 2 to analyze every other record
-delrec 1
-sdate - Start Date for data to be analyzed.
-sdate 01jan70
-stime - Start time (hours, minutes, seconds) for data to be
analyzed. stime=(0,0,0) is 0000 UT on sdate.
-stime 0 0 0
-etime - END time (hours, minutes, seconds) for data to be
analyzed. etime=(0,0,0) is 0000 UT on sdate.
-etime 999 59 59
-ctrlfile - Inscal control file. This contains parameters which
control the ACF fitting process.
-ctrlfile insc_ctrl
-parmsfile - Inscal calibration file.
-parmsfile mhr_parms.dat
Example: /usr/oasis/oasis2004/Oasis/inscal/source/inscal -cmdfile insc_cmd_rt_50kHz
In this example all required arguments are in insc_cmd_rt_50kHz and the only inscal argument
specifies the argument file.
insc_cmd_rt_50kHz:
-cmdfile insc_cmd_rt_50kHz
-oasisdir /usr/oasis/experiments/2004/06dec04/
-tag 50-2
-outfile #dummy
-insfile /usr/oasis/experiments/midasw_realtime/status/current_50kHz_cmst.status
-infile #dummy
-runmode realTime
-startrec 1577
-endrec 99999999
-delrec 1
-sdate 06dec04
-stime 0 0 0
-etime 999 59 59
-ctrlfile /usr/oasis/oasis2004/Oasis/inscal/controlFiles/insc_ctrl_diag
-parmsfile /usr/oasis/oasis2004/Oasis/inscal/controlFiles/mhr_parms.dat_0
VER - INSCAL version number (*100). VER=800 for version 8.0.
FMT='(A3, 3X, I4)'
QSSW - Control switches
FMT='(A4, 3X, 16(I3))'
QSSW( 1) -
QSSW( 2) - QOUT1=1 - PRINT LEVEL 1 DIAGNOSTIC OUTPUT
(LEAST DETAILED)
QSSW( 3) - QOUT2=1 - PRINT LEVEL 2 DIAGNOSTIC OUTPUT
(INTERMEDIATE DETAIL)
QSSW( 4) - QOUT3=1 - PRINT LEVEL 3 DIAGNOSTIC OUTPUT
(MOST DETAILED)
QSSW( 5) - QOUTT=1 - PRINT TERMINAL DIAGNOSTIC OUTPUT
QSSW( 6) - QFTZL=1 - INCLUDE ZERO LAG PRODUCT IN FIT
QSSW( 7) - QFFIT=1 - SEARCH ALL ACFS REGARDLESS OF SNR
QSSW( 8) - QQUIK=1 - USE QUICK SEARCH FOR INITIAL GUESS
QSSW( 9) - QGRID=1 - USE GRID SEARCH FOR INITIAL GUESS
QSSW(10) - QRCOR=1 - APPLY RANGE CORRECTION
USING ASET(11)
QSSW(11) - QACOR=1 - APPLY SPECTRAL ASSYMETRY CORRECTION
USING ASET(8)
QSSW(12) - QRTIM=1 - FOLLOW AND FIT REAL-TIME DATA FILE
QSSW(13) - QISDB=1 - UPDATE INSCOHERENT SCATTER DATABASE
WITH RESULTS OF FIT
QSSW(14) - QERFP=1 - DO E-REGION FULL PROFILE FIT
QSSW(15) - QRCPL=1 - MAKE PLOTS OF FIT PARAMETERS
QSSW(16) - QNPTF>0 - MULTIPLY NO. OF POINTS IN SPECTUM CALCULATION BY QNPTF
QSSW(17) - QRTUT=1 - MAKE RADAR TUTORIAL
QSSW(18) - QFTFN=1 - FIT TI, TR
2 - FIT (TE+TI)/2, (TE-TI)/2
3 - FIT TI, TE
QSSW(19) - QIFMI=1 - INPUT DATA WITH IFAM_IN
QSSW(20) - QIFMO=1 - OUTPUT DATA WITH IFAM_OUT
QSSW(21) -
QSSW(22) -
QSSW(23) -
QSSW(24) -
NMODEA - Number of modes to analyze. If NMODES<0, all modes are
analyzed. If NMODES>=0, only the modes specified by MODE
are analyzed. Up to 32 modes may be specified.
FMT='(A6, 3X, I2)'
MODEA - List of modes to be analyzed. These are compared with
the MODE parameter in ifam (see ifamf.h).
FMT='(A32)'
PARMF - Fixed parameters. Up to 16 may be specified.
FMT='(E13.5)'
PARMF( 1) = Calculate theoretical ACF to at least
this lag
NZR - Number of altitude segments to be analyzed. Search
parameters are different in each altitude segment.
Ther may be at most 6 altitude segments.
FMT='(A3, 3X, I1)'
Z1 - Lowest altitude in an altitude segment.
FMT='(A2, 3X, 6(F6.1,1X))'
Z2 - Highest altitude in an altitude segment.
FMT='(A2, 3X, 6(F6.1,1X))'
SMLIM - Lower limit on signal to noise ratio. Data for which
SNR
Profile integration specification file
For a specified mode, include a file of the form modename.dat in the directory (oasisdir) in which inscal
is run, e.g. a480.30m8910.dat. This file has the following format:
MDNAME = A480.30M8910
MDPTYPE = 2
MDNROWS = 5
MDNCOLS = 3
150.0 200.0 3.0
200.0 300.0 6.0
300.0 400.0 8.0
400.0 500.0 10.0
500.0 800.0 12.0
MDPTYPE = 1 (N-point centered mean), MDTYPE = 2 (N-point
centered median), MDTYPE = 3 (L2 spline fit) AMD MDPTYPE= 4 (L1
spline fit) are supported. the integration parameters are
mdptype dependent and are entered as an array with mdrows rows
and mdcols columns.
For MDPTYPE = 1 and 2, each row contains a minimum and maximum
altitude and the a specification of the number of points in the
mean or median. thus the first row states that between 200 and
300 km a 13-point average will be computed. acfs below 150 km
are left unchanged in this example.
For MDTYPE = 3 and 4, the format is as follows
MDNAME = A480.30M8910
MDPTYPE = 4
MDNROWS = 1
MDNCOLS = 4
100.0 5.0 600.0 100.0
In the example,
Z1 = 100.0 - Lowest spline knot
DZ1 = 5.0 - Distance from z1 to next spline knot
Z2 = 600.0 - Highest spline knot
DZ2 = 100.0 - Approximate knot spacing at z2
In practice the actual spacing between the second-highest knot
and z2 will be somewhat larger than dz2. Knot spacing increases
linearly with altitude. the example yields the following knots:
100.0 105.0 111.2 118.6 127.5 138.1 150.9 166.2 184.5
206.4 232.6 264.0 301.7 346.8 400.8 465.5 600.0
Verifying Inscal
A simple regression analysis is performed by "make verify". Results
are in verify.txt. Changes in inscal or its input parameters may
change the output files prduced by make verify. Also, the inscal list
file (insc_list_v for these checks) contains several lines which always
change. So, verify inscal is an aid to verifying changes in inscal, but
is not an automatic validity checker. The output files in the CVS
archive, insc_list_v_rock, mlh_v_rock and print_mlh_v_rock may
need to be updated, following thorough checking, after inscal has been
modified.
Running inscal in real-time mode
Inscal can analyze data in real-time. In this case it typically is executed from
within a script.
Inscal control input is organized so that all parameters which typically need to be
changed from experiment to experiment are command line parameters. It is also possible
to read these parameters from a command file. Basically there are two ways of running
real-time inscal
1. Read all parameters from a file.
The following example shows how the 28 January, 2005 realtime analyis
was run in this manner for single-pulse data:
/usr/oasis/oasis2004/Oasis/inscal/bin/inscal -cmdfile insc_cmd_rt_50kHz_28jan05
where insc_cmd_rt_50kHz_28jan05 was:
-cmdfile insc_cmd_rt_50kHz_28jan05
-oasisdir /usr/oasis/experiments/2005/28jan05/
-tag #startTime
-outfile #dummy
-insfile /usr/oasis/experiments/midasw_realtime/status/current_50kHz_cmst.status
-infile #dummy
-runmode realTime
-startrec 1
-endrec 99999999
-delrec 1
-sdate 28jan05
-stime 0 0 0
-etime 999 59 59
-ctrlfile /usr/oasis/oasis2004/Oasis/inscal/controlFiles/insc_ctrl_0
-parmsfile /usr/oasis/oasis2004/Oasis/inscal/controlFiles/mhr_parms.dat_0
To run another day, change oasisdir, insfile and sdate as appropriate.
When read from a command file, cmdfile is ignored. Inscal writes a
command file containing the commands actually used to oasisdir,
including the value of cmdfile.
2. Read default parameters from a file and overide some of these on the command line.
For example:
/usr/oasis/oasis2004/Oasis/inscal/bin/inscal -cmdfile /usr/oasis/oasis2004/Oasis/inscal/controlFiles/insc_cmd_rt -oasidir /usr/oasis/experiments/2005/28jan05/ -insfile
/usr/oasis/experiments/midasw_realtime/status/current_50kHz_cmst.status -sdate 28jan05
When runmode=realTime, inscal gets the name of the output Madrigal file
from insfile. This will usually be a Madrigal experiment directory. In
addition to the Madrigal file, inscal places a link to inscal output
files and record plots in the Madrigal experiment directory. The
inscal output files are written to oasisdir, which should normally be
of the form shown in the example. Inscal creates oasisdir and the
Madrigal experiment directory and its subdirectories if they do not
already exist. However, the experiment will not be accessible through
Madrigal until it has been properly entered into Madrigal.
When running in real-time mode, inscal will continue polling insfile
indefinitely. For most experiments, inscal should be stopped at the
conclusion of the experiment. SIGTERM will signal inscal to exit
gracefully.
After analyzing each record, inscal writes a one-line status file of
the form inscalStatusFile_tag to oasisdir. The fields are
defined by the following Fortran write statement:
UTIME = TIME()
PID = GETPID()
WRITE (LFN,FMT=
*'(I6, 2X, A, 2X, I10, 2X, I8, 2X,
*I2.2, 1X, I2.2, 1X, I4.4, 2X,
*I2.2, 1X, I2.2, 1X, I2.2, 2X,
*I2.2, 1X, I2.2, 1X, I2.2,
*I4, 4F7.1, 2X, A)') PID,OUTFILE(1:LNBLNK(OUTFILE)),UTIME,RECNUM,
* IDATE(1,ICACF),IDATE(2,ICACF),IDATE(3,ICACF),
* STIME(1,ICACF),STIME(2,ICACF),STIME(3,ICACF),
* ETIME(1,ICACF),ETIME(2,ICACF),ETIME(3,ICACF),
* KINST,
* BAZ(ICACF),EAZ(ICACF),BEL(ICACF),EEL(ICACF),
* MODE(1) (1:LNBLNK(MODE(1)))
>
where the output parameters are from ifam, except for utime, pid and
recnum. Recnum is as defined in the Midas_io Fortran Library, i.e.
recnum=0 for the first record.
For example:
4786 /opt/madrigal/experiments/2005/mlh/28jan05/mlh050128a.000 1107441325 8 28 01 2005 16 14 00 16 18 00 32 178.0 178.0 88.0 88.0 s480.20z8910
This file can be polled to determine whether inscal is running
properly. In the examples, inscal used its Unix start time as the tag.
Altenatively. a controlling script might specify tag to inscal, for
example, the Unix time when it initiated inscal. This would ensure that
the script knows the name of the appropriate inscal status file.
analyzeData
THE FOLLOWING IS OUT OF DATE
When not in real-time mode, typically Inscal is run by the analyzeData Tcl script.
Usage: analyzeData [options]
Options
-inscal inscal_executable
-expStartDate experiment_start_date
-inscalCtrl inscal_control_file
-mhrParms mhr_parms.dat_template
-calibrate calibration_executable
-calibrateCtrl calibration_control_file_template
-expDir experiment_directory
-expLabel experiment_label
-doCal do_calibration_flag
-dataFile data_file
Example:
/usr/oasis/inscal2002/Oasis/inscal/bin/analyzeData -expDir /usr/oasis/inscal2002/midaswCalibration/20feb04 -dataFile /usr/oasis/experiments/2004/20feb04/run/2004-02-20_50kHz.ifam
If no data files are specified, this help message is printed.
If one or more data files are specified, they are analyzed by inscal.
Data records will be anlyzed by inscal if they begin after incsal control file parameter
STIMEC (UT) relative to expStartDate and end before ETIMEC (UT) relative to
experimentStartDat.
If calibration is specified by the doCal option, an electron density calibration
is calculated.
-inscal - inscal I.S. analysis program executable
default = /usr/oasis/inscal2002/Oasis/inscal/bin/inscal
-expStartDate - experiment start date, e.g. 23jun05
default = 01jan02
-inscalCtrl - inscal control file template
default = /usr/oasis/inscal2002/Oasis/inscal/controlFiles/insc_ctrl_0
-mhrParms - mode parameter file for inscal
default = /usr/oasis/inscal2002/Oasis/inscal/controlFiles/mhr_parms.dat_0
-calibrate - calibration program executable
default = /usr/oasis/inscal2002/Oasis/inscal/bin/calibrate
-calibrateCtrl - calibration program control file
default = /usr/oasis/inscal2002/Oasis/inscal/controlFiles/calibrate.cmd
-expDir - experiment directory for results; it will be created if necessary
default = inscal_UnixTimeInSeconds
-expLabel - string used in plot titles and labels
default = Inscal run at UnixTimeInSeconds
-doCal - 1 for calibration experiment, otherwise, 0
default = 0
-dataFile - the data files: may be more than one; each will be analyzed in turn
default = None
INSCAL Call Graph
1 INSCAL
2 GETWTIME [external]
3 CEDAROPEN [external]
4 GETCWD [external]
5 HOSTNM [external]
6 LNBLNK [external]
7 OPENMHRPARMS [external]
8 OPENMIDASFILE [external]
9 READIFAMRECORD [external]
10 CKIFAM
11 ERROR
12 FLUSH [external]
13 GETCMDFILE
14 GETCTRLFILE
15 GETDELREC
16 GETENDREC
17 GETETIME
18 LNBLNK [external]
19 GETINFILE
20 GETINSFILE
21 GETOASISDIR
22 GETOUTFILE
23 GETPARMSFILE
24 GETRUNMODE
25 GETSDATE
26 GETSTARTREC
27 GETSTIME
28 LNBLNK [external]
29 GETTAG
30 ACFIT
31 FCHISQ
32 LNBLNK [external]
33 AQUIK
34 BNEWT
35 DFUN
36 VELX
37 BADAQ
38 DOPP
39 FITPH
40 CURFIT
41 FCHISQ
42 FDERIV
43 FUNCTLWE
44 ACFMINI
45 ACFP2ACFM
46 ACFPINI
47 ATM2ACFP
48 DFUN
49 FOURS
50 FOUR2
51 BITRV
52 COOL2
53 FIXRL
54 SPECT2
55 DAW
56 SPECTC
57 WOFZ
58 ATMLWE
59 DENSTY
60 ATMPRN
61 FUNCTLWE [see line 43]
62 FUNCTN
63 DFUN
64 FOURS [see line 49]
65 FULLSPEC
66 FOUR1
67 SP2TITE
68 SPECT2 [see line 54]
69 SPECTC [see line 56]
70 FUNPH
71 MINV
72 FUNCTLWE [see line 43]
73 FUNCTLWE [see line 43]
74 FUNCTN [see line 62]
75 FUNPH
76 FUNPH
77 ARGET
78 ARUPD
79 BADACF
80 CONVRT [external]
81 EFITFP
82 DENSTY
83 ACFDPRN
84 ACFMPRN
85 ACFPPRN
86 ATM2OFAM
87 ATMINI
88 ATMMOD
89 ATMMSIS
90 GTS5 [external]
91 METERS [external]
92 TRETRV [external]
93 TSELEC [external]
94 ATMPRN
95 CURFIT [see line 40]
96 FUNCTLWE [see line 43]
97 GTS5 [external]
98 METERS [external]
99 OFAM2ATM
100 SETCON
101 ERROR [see line 11]
102 TRETRV [external]
103 TSELEC [external]
104 EXINSC
105 CLOSEMIDASFILE [external]
106 CEDARCLOSE [external]
107 EXIT [external]
108 FDATE [external]
109 FLUSH [external]
110 IEEE_FLAGS [external]
111 PRPARM
112 LNBLNK [external]
113 GETPARMSFILE
114 PRSTAT [external]
115 FIT1
116 A2SP
117 CURFIT [see line 40]
118 DSP2DA
119 DSP2DT
120 SP2TITE
121 TITE2A
122 FITPH [see line 39]
123 FITSTATS
124 FLUSH [external]
125 IF2OFI
126 PARAM
127 PARMOD
128 ERROR [see line 11]
129 POINT [external]
130 SETCON [see line 100]
131 SETFIT
132 SETSTAT
133 STACOR
134 CONVRT [external]
135 ARGIN
136 GETCWD [external]
137 IARGC [external]
138 LNBLNK [external]
139 TIME [external]
140 EXINSC [see line 104]
141 GETARG [external]
142 READARGS
143 LNBLNK [external]
144 TIME [external]
145 GETOKEN
146 ARINI
147 CEDARCLOSE [external]
148 CTRLIN
149 CTRLOUT
150 EXINSC [see line 104]
151 FDATE [external]
152 FIXIFAM
153 FIXSPEC
154 FOUR1
155 GTCMTM
156 DATER [external]
157 JDAY [external]
158 GETETIME [see line 17]
159 GETSDATE
160 GETSTIME [see line 27]
161 FLUSH [external]
162 GEOPRM
163 GETDATAFILE [external]
164 GETMADFILE [external]
165 GETMIDASIOERR [external]
166 GTCMTM [see line 155]
167 GTIFTM
168 JDAY [external]
169 IFAMINF
170 IFAMOUTNF
171 INICOM
172 MAKENOTES
173 LNBLNK [external]
174 GETOASISDIR
175 GETETIME [see line 17]
176 GETOUTFILE
177 GETSDATE
178 GETSTIME [see line 27]
179 GETTAG
180 SPLITMF
181 LNBLNK [external]
182 OFAM2CEDAR
183 CEDARGETSCALEFACTOR [external]
184 CEDARCREATE [external]
185 CEDARSET1DPARM [external]
186 CEDARSET2DPARM [external]
187 CEDARWRITE [external]
188 JDATER [external]
189 JDAY [external]
190 PLOTPAR [external]
191 POSTINT
192 SHELL1
193 PRECSM
194 LNBLNK [external]
195 FLUSH [external]
196 PROFINT
197 DEFINELPFIT
198 KNOTSION
199 FTL1L2
200 COVC
201 COV
202 INCLDC
203 SROTM
204 SROTMG
205 L1
206 PSLV
207 REGRSC
208 SHELL1
209 LPBASIS
210 BSPLVD
211 BSPLVN
212 INTERV
213 GTMDAT
214 KNOTSION
215 SHELL1
216 LNBLNK [external]
217 PRSFIL
218 GETPID [external]
219 LNBLNK [external]
220 TIME [external]
221 FLUSH [external]
222 SETOUTFILE
223 SETSIGNALS [external]
224 SPLITMF [see line 180]
225 SYSTEM [external]
226 WRITEARGS
227 LNBLNK [external]
228 EXINSC [see line 104]
229 SYSTEM [external]
230 CEDAR2OFAM
231 GRECNO [external]
232 MDGKST [external]
233 MDGNRW [external]
234 MDIBCS [external]
235 MDIBDT [external]
236 MDIBHM [external]
237 MDIBYR [external]
238 MDIECS [external]
239 MDIEDT [external]
240 MDIEHM [external]
241 MDIEYR [external]
242 MDISDR [external]
243 JDAY [external]
244 MDG1DP [external]
245 MDGEOD [external]
246 MDGFLT [external]
247 GETEXT
248 LNBLNK [external]
249 PRNTEXT
250 LNBLNK [external]
251 IFAMOUTF
252 POLBAS
253 LEGBAS
254 HRMBAS
255 SPLBASI
256 BSPLVD [see line 210]
257 INTERV
258 SPLBAS22
259 SPLBAS
260 BSPLVD [see line 210]
261 INTERV
262 SPLBASP2
263 SPLBAS [see line 259]
264 SPLBASP
265 BSPLVD [see line 210]
266 INTERV
267 KNOTS1
268 KNOTS2
269 KNOTSP
270 PARINI
271 SP2A
272 TITE2SP
273 A2TITE
INSCAL Procedure Synopsis
PROGRAM INSCAL ()
C
C INCOHERENT SCATTER DATA ANALYSIS PROGRAM
C JOHN M. HOLT
C MIT HAYSTACK OBSERVATORY
C WESTFORD, MASSACHUSETTS 01886
C
C VERSION 1.0 - 10/83 - ORIGINAL HARRIS VERSION
C VERSION 2-6 - SEE HARRIS DOCUMENTATION
C VERSION 6.0 - 12/25/92 - MIDAS/IFAM/OFAM/MADRIGAL PRODUCTION
C VERSION 6.1 - 07/16/93 - NUMEROUS SMALL ENHANCEMENTS
C VERSION 7.0 - 11/22/93 - FULL-PROFILE E-REGION SEARCHES
C VERSION 7.1 - 06/06/94 - OUTPUTS HTML FILE FOR MOSAIC
C VERSION 8.0 - 11/21/97 - SOLARIS VERSION W ASSYMETRY CORRECTION
C VERSION 8.1 - 10/04/99 - DOPPLER VELOCITY, STATISTICS
C VERSION 8.2 - 08/22/00 - MADREC OUTPUT, RT SUPPORT
C VERSION 9.0 - 03/07/03 - MIDASW SUPPORT, MADLIB FREE
C VERSION 10.0 - 06/??/04 - FULL MIDASW SUPPORT, MULTIPLE ACFS
C
C INSCAL ANALYZES INCOHERENT SCATTER AUTOCORRELATION FUNCTIONS
C (ACFS) TO DETERMINE IONOSPHERIC PLASMA PARAMETERS. THE INPUT
C DATA ARE TAKEN FROM IFAMF. THE OUTPUT PARAMETERS ARE PLACED
C IN OFAMF. INSCAL CONTROL PARAMETERS ARE IN CNTRL. DATA FROM
C ANY I. S. RADAR MAY BE ANALYZED BY WRITING A ROUTINE TO CONVERT
C THE NATIVE DATA FORMAT OF THAT RADAR TO IFAMF FORMAT. OUTPUT
C PARAMETERS MAY BE SAVED IN ANY FORMAT BY WRITING A ROUTINE TO
C CONVERT OFAMF TO THE DESIRED OUTPUT FORMAT.
C
C IN THIS VERSION OF INSCAL, THE NATIVE RAW DATA FORMAT IS THE
C MIDAS CMST STRUCTURE USED AT MILLSTONE HILL, AND THE OUTPUT
C FORMAT IS THE NCAR DATABASE MADRIGAL FORMAT. IN ADDITION, THE
C NCAR GRAPHICS PACKAGE IS USED TO PRODUCE A SUMMARY PLOT OF EACH
C ANALYZED RECORD.
C
C Use egrep '(\.\.\.\.\.)' inscal.f for INSCAL summary.
C
C Usage: inscal control_file working_OASISDIR
C [file_tag output_file real_time_status_file inscal_status_file]
C
SUBROUTINE ACFDPRN (LFNPRN,ALT1,ALT2,FAC)
INTEGER LFNPRN !Read-only
DOUBLE PRECISION ALT1 !Read-only
DOUBLE PRECISION ALT2 !Read-only
DOUBLE PRECISION FAC(*) !Read-only
C
C ACFDPRN PRINTS THE DIFFERENCE BETWEEN THE MEASURED
C AUTOCORRELATION FUNCTIONS IF IFAM AND ACFM.
C
C LFNPRN - LFN FOR OUTPUT
C ALT1 - BEGIN ALTITUDE TO PRINT
C ALT2 - END ALTITUDE TO PRINT
C FAC - FACTOR TO APPLY TO ACFM ACF BEFORE TAKING DIFFERENCE
C
C JMH - 11/97
C
SUBROUTINE ACFIT (NZR,Z1,Z2,NSI,SNLIM,NITER,ASET,IER)
INTEGER NZR !Read-only
DOUBLE PRECISION Z1(6) !Read-only
DOUBLE PRECISION Z2(6) !Read-only
INTEGER NSI(NFPMAX,6) !Read-only
DOUBLE PRECISION SNLIM(6) !Read-only
INTEGER NITER(6) !Read-only
DOUBLE PRECISION ASET(NFPMAX,6) !Read-only
INTEGER IER
C
C JMH - 11/83 HARRIS FORTRAN 77
C JMH - 12/92 Initial production Sun Version
C
C ACFIT FITS THEORETICAL ACFS TO THE EXPERIMENTAL ACFS PASSED TO IT
C IN COMMON/IFAM/. IF THE SEARCH IS SUCCESSFUL, THE VALUES PUT IN A
C BY THE SUCCESSFUL FIT ARE UNORMALIZED AND PUT IN AA. FINALLY,
C THE VALUES IN AA, WHETHER "BAD" OR "GOOD", ARE TRANSFERED TO PART
C AND HENCE TO THE DATABASE RECORD.
C Use egrep '(\.\.\.\.\.)' acfit.f for ACFIT summary
C
SUBROUTINE ACFMINI ()
C
C ACFMINI INITIALIZES THE TAU ARRAY IN ACFM
C
C JMH - 11/93
C
SUBROUTINE ACFMPRN (LFNPRN,FAC)
INTEGER LFNPRN !Read-only
DOUBLE PRECISION FAC(*) !Read-only
C
C ACFMPRN PRINTS THE MEASURED PART OF THE ACF IN ACFM
C
C JMH - 11/93
C
SUBROUTINE ACFP2ACFM ()
C
C JMH - 11/93 ADAPTED FROM FUNCTN
C
C FUNCTN CALCULATES AN AUTOCORRELATION OBJECT FUNCTION MATRIX FROM
C AN ATMOSPHERIC PARAMETERS MATRIX.
C
SUBROUTINE ACFPINI ()
C
C JMH - 11/93
C
SUBROUTINE ACFPPRN (LFNPRN)
INTEGER LFNPRN !Read-only
C
C JMH - 11/93
C
SUBROUTINE AQUIK (NLAG,DTAU,PO,ACFR,ACFI,NE,TI,TR,VO,IER)
INTEGER NLAG !Read-only
DOUBLE PRECISION DTAU !Read-only
DOUBLE PRECISION PO !Read-only
DOUBLE PRECISION ACFR(*) !Read-only
DOUBLE PRECISION ACFI(*) !Read-only
DOUBLE PRECISION NE
DOUBLE PRECISION TI
DOUBLE PRECISION TR
DOUBLE PRECISION VO
INTEGER IER
C
C JMH - 11/82 HARRIS FORTRAN 77
C
C AQUIK COMPUTES ION AND ELECTRON TEMPERATURE FROM THE FIRST ZERO
C CROSSING AND DEPTH OF FIRST MINIMUM OF AN ACF, AND THE VELOCITY
C FROM THE PHASE OF THE ACF. THYPICALLY AQUICK PROVIDES A FIRST
C GUESS OF THE PARAMETERS PRIOR TO A NON-LINEAR LEAST-SQUARES FIT.
C
C INPUT PARAMETERS
C NLAG - NUMBER OF ACF LAGS
C DTAU - ACF LAG SPACING (MICRO SEC)
C PO - LOG POWER, NORMALIZED TO NE=.5*(1+TR+DF)*(1+DF)*10**PO
C ACFR - REAL PART OF ACF, NORMALIZED TO 1.0 AT ZERO LAG
C ACFI - IMAGINARY PART OF ACF
C
C OUTPUT PARAMETERS
C NE - ELECTRON DENSITY (1/CM**2)
C TI - ION TEMPERATURE (DEG K)
C TR - ELECTRON/ION TEMPERATURE RATIO
C VO - LINE-OF SIGHT DRIFT (M/SEC)
C IER=0 - SUCCESSFUL NE,TI,TR,VO DETERMINATION
C 1 - FIRST TWO LAGS NOT BOTH > 0
C 2 - ALL ACFR POINTS > 0
C 3 - IMAGINARY ZERO CROSSING
C 4 - NO MINIMUM FOUND AFTER FIRST ZERO CROSSING
C 5 - FIRST MINIMUM OUT OF TABULATED RANGE
C 6 - SUCCESSFUL TI,TR,VO DETERMINATION, PO OUT OF RANGE
C 7 - TI OR TR OUT OF RANGE
C 8 - VO OUT OF RANGE
C
SUBROUTINE ARGIN (ISTAT)
INTEGER ISTAT
C
C JMH - 10/25/04
C
C ARGIN PARSES INSCAL COMMAND LINE ARGUMENTS.
C
SUBROUTINE READARGS (ISTAT)
INTEGER ISTAT
C
C JMH - 11/05/04
C
C READARGS READS INSCAL COMMAND LINE ARGUMENTS FROM A FILE
C
SUBROUTINE WRITEARGS (CMDFILEOUT,ISTAT)
CHARACTER CMDFILEOUT !Read-only
INTEGER ISTAT !Read-only
C
C JMH - 11/11/04
C
C WRITEARGS WRITES INSCAL COMMAND LINE ARGUMENTS TO A FILE
C
SUBROUTINE SETOUTFILE (OUTFILEIN)
CHARACTER * 128 OUTFILEIN !Read-only
C
SUBROUTINE ARINI ()
C
SUBROUTINE ARUPD (Z,TI,DTI,TR,DTR)
DOUBLE PRECISION Z !Read-only
DOUBLE PRECISION TI !Read-only
DOUBLE PRECISION DTI !Read-only
DOUBLE PRECISION TR !Read-only
DOUBLE PRECISION DTR !Read-only
C
SUBROUTINE ARGET (Z,TI,TR,IER)
DOUBLE PRECISION Z !Read-only
DOUBLE PRECISION TI
DOUBLE PRECISION TR
INTEGER IER
C
SUBROUTINE ATM2ACFP ()
C
C JMH - 11/93 ADAPTED FROM FUNCTN
C
C ATM2ACFP CALCULATES AN AUTOCORRELATION OBJECT FUNCTION MATRIX FROM
C AN ATMOSPHERIC PARAMETERS MATRIX.
C
SUBROUTINE ATM2OFAM ()
C
C JMH - 06/94
C
SUBROUTINE ATMINI ()
C
C JMH - 10/93
C
SUBROUTINE ATMLWE (ZL,ZU,DLB,TINF,TLB,XM,ZLB,S2,T0,ZA,Z0,TR12,CLF,CPO,CVO)
DOUBLE PRECISION ZL !Read-only
DOUBLE PRECISION ZU !Read-only
DOUBLE PRECISION DLB !Read-only
DOUBLE PRECISION TINF !Read-only
DOUBLE PRECISION TLB !Read-only
DOUBLE PRECISION XM !Read-only
DOUBLE PRECISION ZLB !Read-only
DOUBLE PRECISION S2 !Read-only
DOUBLE PRECISION T0 !Read-only
DOUBLE PRECISION ZA !Read-only
DOUBLE PRECISION Z0 !Read-only
DOUBLE PRECISION TR12 !Read-only
DOUBLE PRECISION CLF !Read-only
DOUBLE PRECISION CPO(*) !Read-only
DOUBLE PRECISION CVO(*) !Read-only
C
C Z1 - lower altitude boundry for calculation altitude (km)
C Z2 - upper altitude boundry for calculation altitude (km)
C DLB - density at lower boundry ZLB
C TINF - exospheric temperature (K)
C TLB - temperature at ZLB (K)
C XM - molecular weight of gas species
C ZLB - lower boundry (120 km)
C S2 - TLB'/(TINF-TLB)
C T0 - mesopause temperature (K)
C ZA - altitude of temperature profile junction
C Z0 - mesopause height (km)
C TR12 - mesopause shape parameter
C CLF - log collision frequency factor
C CPO - polynomial coefficients for PO
C
SUBROUTINE ATMMOD ()
C
C JMH - 10/93
C
SUBROUTINE ATMMSIS (ZL,ZU)
DOUBLE PRECISION ZL !Read-only
DOUBLE PRECISION ZU !Read-only
C
SUBROUTINE ATMPRN (LFNPRN,INDPRN)
INTEGER LFNPRN !Read-only
INTEGER INDPRN(*) !Read-only
C
C JMH - 10/93
C
SUBROUTINE BADACF (A,SIGA,ISTBAD)
DOUBLE PRECISION A(*)
DOUBLE PRECISION SIGA(*)
INTEGER ISTBAD
C
C JMH - 9/81
C
C BADACF SETS ACF SEARCH PARAMETERS TO A STANDARD BAD VALUE.
C
C A - SEARCH PARAMETERS
C SIGA - SEARCH PARAMETER ERRORS
C ISTBAD - ISTBAD=1 -> BAD PARAMETERS
C
SUBROUTINE BADAQ (TI,TR,VO)
DOUBLE PRECISION TI
DOUBLE PRECISION TR
DOUBLE PRECISION VO
C
C JMH - 11/82
C
DOUBLE PRECISION FUNCTION BNEWT (I,X,DX,Y,XVAL)
INTEGER I !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION DX !Read-only
DOUBLE PRECISION Y(*) !Read-only
DOUBLE PRECISION XVAL !Read-only
C
C JMH - 11/82
C
SUBROUTINE CEDAR2OFAM (MADFIL,RECNO,IER)
INTEGER MADFIL !Unknown (program incomplete)
INTEGER RECNO !Unknown (program incomplete)
INTEGER IER
C
C CEDAR2OFAM reads a cedar record from an in-memory madrec structure
C and populates an ofam structure.
C
C Inputs: INTEGER MADFIL - pointer to in-memory madrec structure
C (that is, opened with iotype 30 or greater)
C
C INTEGER RECNO - record number of interest
C
C Effect: Populates ofam.h common blocks OFAM, OFAMC
C Sets IER to 0 if success, negitive otherwise
C Sets current record in MADFIL to recno
C
C Errors: If MADFIL is NULL, or RECNO does not point to data record
C
LOGICAL FUNCTION CKIFAM (IREC)
INTEGER IREC !Read-only
C
C JMH - 12/92
C
SUBROUTINE CTRLIN (CTRLFIL,ISTAT)
CHARACTER CTRLFIL !Read-only
INTEGER ISTAT
C
C JMH - 1/5/93
C
SUBROUTINE CTRLOUT (CTRLFIL,ISTAT)
CHARACTER CTRLFIL !Read-only
INTEGER ISTAT
C
C JMH - 12/98 - UPDATED TO INSCAL 8.0
C
SUBROUTINE CURFIT (X,Y,NPTS,NTERMS,W,A,SIGMA,FLAMDA,YFIT,CHISQR,CORMAT,NTS,
COVMAT,DERIV,INIT,NS,FUNCTN,NTYPE,MODE)
DOUBLE PRECISION X(*)
DOUBLE PRECISION Y(*) !Read-only
INTEGER NPTS !Read-only
INTEGER NTERMS !Read-only
DOUBLE PRECISION W(*) !Read-only
DOUBLE PRECISION A(*)
DOUBLE PRECISION SIGMA(*)
DOUBLE PRECISION FLAMDA
DOUBLE PRECISION YFIT(*)
DOUBLE PRECISION CHISQR
DOUBLE PRECISION CORMAT(*)
INTEGER NTS
DOUBLE PRECISION COVMAT(*)
DOUBLE PRECISION DERIV(NPTS,*)
INTEGER INIT
INTEGER NS(NTERMS) !Read-only
SUBROUTINE FUNCTN
Argument 1: INTEGER
Argument 2: DOUBLE PRECISION ARRAY
Argument 3: INTEGER ARRAY
Argument 4: INTEGER
Argument 5: DOUBLE PRECISION ARRAY
Argument 6: INTEGER ARRAY
Argument 7: DOUBLE PRECISION ARRAY
INTEGER NTYPE(NPTS) !Read-only
INTEGER MODE !Read-only
C
C SUBROUTINE CURFIT
C
C PURPOSE
C MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION
C WITH A LINEARIZATION OF THE FITTING FUNCTION
C
C USAGE
C CALL CURFIT(X,Y,NPTS,NTERMS,W,A,SIGMA,FLAMDA,YFIT,CHISQR,
C * CORMAT,NTS,COVMAT,DERIV,INIT,NS,FUNCTN,NTYPE,MODE)
C
C INPUT PARAMETERS
C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
C NPTS - NUMBER OF PAIRS OF DATA POINTS
C NTERMS - NUMBER OF SEARCH PARAMETERS
C W - ARRAY OF WEIGHTS FOR Y DATA POINTS
C FLAMDA - PROPORTION OF GRADIENT SEARCH INCLUDED
C NS - NS(I)=1 -> SEARCH PARAMETER A(I)
C FUNCTN - LEAST SQUARES OBJECTIVE FUNCTION
C NTYPE - PASSED THROUGH TO FUNCTN
C MODE - MODE=2 -> DIAGNOSTIC OUTPUT
C INIT - ITERATION COUNTER
C OUTPUT PARAMETERS
C SIGMA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A
C YFIT - ARRAY OF CALCULATED VALUES OF Y
C CHISQR - REDUCED CHI-SQUARE FOR FIT
C NTS - NUMBER OF SEARCHED PARAMETERS
C COVMAT - COVARIANCE MATRIX
C CORMAT - CORRELATION MATRIX
C DERIV - DERIVATIVES OF THE FITTING FUNCTION
C INPUT/OUTPUT PARAMETERS
C A - ARRAY OF SEARCH PARAMETERS
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C FUNCTN (X, NPTS, A, YFIT)
C EVALUATES THE FITTING FUNCTION
C FCHISQ (Y, NPTS, NFREE, WEIGHT, YFIT)
C EVALUATES REDUCED CHI-SQUARE FOR FIT TO DATA
C FDERIV (X, NPTS, A, NTERMS, YFIT, DERIV)
C EVALUATES THE DERIVATIVES OF THE FITTING FUNCTION
C MINV (ARRAY, NTERMS, DET, L, M)
C INVERTS A SYMETRIC TWO-DIMENSIONAL MATRIX
C OF DEGREE NTERMS AND CALCULATED ITS DETERMINANT
C
C COMMENTS
C SET FLAMDA = 0.001 AT BEGINNING OF SEARCH
C
C FROM:
C PHILIP R. BEVINGTON
C DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES
C MCGRAW-HILL BOOK COMPANY, 1969
C
DOUBLE PRECISION FUNCTION DAW (X)
DOUBLE PRECISION X !Read-only
C
C DAW COMPUTES DAWSON'S INTEGRAL, DAW(X), TO 11 SIGNIFICANT
C FIGURES, USING RATIONAL CHEBYSHEV APPROXIMATIONS
C (W. J. CODY ET AL., MATH COMP. 24, 1970, 171-178).
C
SUBROUTINE DEFINELPFIT (X,Y,DY,CB,M,N,TIN)
DOUBLE PRECISION X(MDIM)
DOUBLE PRECISION Y(MDIM)
DOUBLE PRECISION DY(MDIM)
DOUBLE PRECISION CB(MDIM,5)
INTEGER M
INTEGER N
DOUBLE PRECISION TIN(NDIM)
C
C JMH - 1/2004
C
DOUBLE PRECISION FUNCTION DENSTY (ALT,DLB,TINF,TLB,XM,ALPHA,TZ,ZLB,S2,T0,
ZA,Z0,TR12)
DOUBLE PRECISION ALT !Read-only
DOUBLE PRECISION DLB !Read-only
DOUBLE PRECISION TINF !Read-only
DOUBLE PRECISION TLB !Read-only
DOUBLE PRECISION XM !Read-only
DOUBLE PRECISION ALPHA !Read-only
DOUBLE PRECISION TZ
DOUBLE PRECISION ZLB !Read-only
DOUBLE PRECISION S2 !Read-only
DOUBLE PRECISION T0 !Read-only
DOUBLE PRECISION ZA !Read-only
DOUBLE PRECISION Z0 !Read-only
DOUBLE PRECISION TR12 !Read-only
C
C Calculate Temperature and Density Profiles for MSIS models
C
C ALT - altitude (km)
C DLB - density at lower boundry ZLB
C TINF - exospheric temperature (K)
C TLB - temperature at ZLB (K)
C XM - molecular weight of gas species
C ALPHA - thermal diffusion coefficient
C TZ - temperature at ALT (K) (OUTPUT PARAMETER)
C ZLB - lower boundry (120 km)
C S2 - TLB'/(TINF-TLB)
C T0 - mesopause temperature (K)
C ZA - altitude of temperature profile junction
C Z0 - mesopause height (km)
C TR12 - mesopause shape parameter
C
DOUBLE PRECISION FUNCTION DFUN (TR,PFAC)
DOUBLE PRECISION TR !Read-only
DOUBLE PRECISION PFAC !Read-only
C
C DFUN CALCULATES THE DEBYE FACTOR CORRESPONDING TO GIVEN
C TEMPERATURE RATIO AND POWER FACTOR PFAC=4761.*TE*RK**2/(PNORM*P).
C NEWTON'S METHOD IS USED TO SOLVE DFUN*(1+TR+DFUN)*(1+DFUN)-PFAC=0.
C
SUBROUTINE DOPP (ACFR,ACFI,DTAU,NLAGS,DLAG,LMULTI,VEL,IERROR)
DOUBLE PRECISION ACFR(*) !Read-only
DOUBLE PRECISION ACFI(*) !Read-only
DOUBLE PRECISION DTAU !Read-only
INTEGER NLAGS !Read-only
INTEGER DLAG !Read-only
LOGICAL LMULTI !Read-only
DOUBLE PRECISION VEL
INTEGER IERROR
C
C Created 8/3/89 ... DT
C
C F:DOPP computes the Doppler velocity from the average lag-to-lag
C phase change, in contrast to F:VELX, which uses the average phase
C change referred to lag 0.
C
C Input parameters:
C
C ACFR(): Real part of a.c.f.
C
C ACFI(): Imaginary part of a.c.f.
C
C DTAU: Lag spacing in usecs
C
C NLAGS: Maximum lag to use in computation. NLAGS is less than
C or equal to the total a.c.f. length
C
C DLAG: Pair separation for phase change. DLAG=1 uses
C (0,1), (1,2), etc, while DLAG=2 uses
C (0,2), (1,3), etc.
C
C LMULTI: Multipulse indicator tells the algorithm to NOT use
C the missing lags 0 and 6.
C
C Output parameters:
C
C VEL: Doppler velocity
C
C IERROR: IERROR = 0 ---> velocity o.k.
C IERROR = 1 ---> error while computing velocity
C
C *** NOTE ***
C
C The constant FACTOR is hard-wired for Millstone Hill, and is
C equal to 10**6 * (wavelength / ( 4 * pi ))
C
SUBROUTINE EFITFP (ALT1,ALT2,IER)
DOUBLE PRECISION ALT1 !Read-only
DOUBLE PRECISION ALT2 !Read-only
INTEGER IER
C
C JMH - 10/93
C
SUBROUTINE ERROR (IER,NAME)
INTEGER IER !Read-only
CHARACTER * 6 NAME !Read-only
C
C ERROR PRINTS AN ERROR MESSAGE WHEN AN ERROR IS ENCOUNTERED BY THE
C CALLING ROUTINE. IER INDICATES THE SEVERITY OF THE ERROR AND IS OF
C THE FORM TYPE+N WHERE
C TYPE = 200 IMPLIES WARNING WITH FIX
C TYPE = 100 IMPLIES WARNING
C N = ERROR CODE DEFINED IN CALLING ROUTINE
C NAME IS A SIX CHARACTER LITERAL STRING CONTAINING THE NAME OF
C THE CALLING ROUTINE.
C
SUBROUTINE EXINSC (SIGNAL)
INTEGER SIGNAL !Read-only
C
C JMH - 3/93
C
C EXINSC EXITS INSCAL GRACEFULLY
C
DOUBLE PRECISION FUNCTION FCHISQ (Y,NPTS,NFREE,W,YFIT)
DOUBLE PRECISION Y(*) !Read-only
INTEGER NPTS !Read-only
INTEGER NFREE !Read-only
DOUBLE PRECISION W(*) !Read-only
DOUBLE PRECISION YFIT(*) !Read-only
C
C PURPOSE
C EVALUATE REDUCED CHI SQUARE FOR FIT TO DATA
C FCHISQ = SUM((Y-YFIT)**2 / SIGMA**2) / NFREE
C
C USAGE
C RESULT = FCHISQ(Y, NPTS, NFREE, W, YFIT)
C
C DESCRIPTION OF PARAMETERS
C Y - ARRAY OF DATA POINTS
C NPTS - NUMBER OF DATA POINTS
C NFREE - NUMBER OF DEGREES OF FREEDOM
C W - ARRAY OF WEIGHTS
C YFIT - ARRAY OF CALCULATED VALUES OF Y
C
C FROM:
C PHILIP R. BEVINGTON
C DATA REDUCTION AND ERROR ANALYSIS FOR THE PHYSICAL SCIENCES
C MCGRAW-HILL BOOK COMPANY, 1969
C
SUBROUTINE FDERIV (NPTS,X,NTYPE,NTERMS,A,NS,YFIT,DERIV,FUNCTN)
INTEGER NPTS !Read-only
DOUBLE PRECISION X(*)
INTEGER NTYPE(*) !Read-only
INTEGER NTERMS !Read-only
DOUBLE PRECISION A(*)
INTEGER NS(*) !Read-only
DOUBLE PRECISION YFIT(*) !Read-only
DOUBLE PRECISION DERIV(NPTS,*)
SUBROUTINE FUNCTN
Argument 1: INTEGER
Argument 2: DOUBLE PRECISION ARRAY
Argument 3: INTEGER ARRAY
Argument 4: INTEGER
Argument 5: DOUBLE PRECISION ARRAY
Argument 6: INTEGER ARRAY
Argument 7: DOUBLE PRECISION Array Element
C
C PURPOSE
C EVALUATE DERIVATIVES OF FUNCTION FOR LEAST-SQUARES SEARCH
C FOR ARBITRARY FUNCTION GIVEN BY FUNCTN
C
C USAGE
C CALL FDERIV(X, NPTS, A, NTERMS, YFIT, DERIV)
C
C DESCRIPTION OF PARAMETERS
C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE
C NPTS - NUMBER OF DATA POINTS
C A - ARRAY OF PARAMETERS
C NTERMS - NUMBER OF PARAMETERS
C YFIT - ARRAY OF FUNCTION VALUES
C DERIV - DERIVATIVES OF FUNCTION
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C FUNCTN (NPTS, X, NTYPE, NTERMS, A, NS, DERIV)
C EVALUATES YFIT
C
SUBROUTINE FOURS (DATA,N,ISIGN,ISYM)
DOUBLE PRECISION DATA(*)
INTEGER N !Read-only
INTEGER ISIGN !Read-only
INTEGER ISYM !Read-only
C FOURIER TRANSFORM OF REAL, SYMMETRIC DATA BY THE FAST FOURIER
C TRANSFORM. WRITTEN BY JAMES COOLEY, IBM, AND NORMAN BRENNER, MIT,
C OCTOBER 1968.
C THE INPUT ARRAY IS REAL, LENGTH N/2+1, REPRESENTING A REAL ARRAY
C OF LENGTH N. N IS AN ARBITRARY MULTIPLE OF FOUR, THO IF IT IS
C NOT A POWER OF TWO, REPLACE FOUR2 BELOW BY FOURT.
C THE FULL ARRAY OF LENGTH N HAS EITHER EVEN OR ODD SYMMETRY (ISYM=
C +1 OR -1) ACCORDING TO WHETHER DATA(N+2-I)=+DATA(I) OR
C - DATA(I), FOR I = 2,3,...N. DATA(1) = + OR - DATA(1). ODD
C SYMMETRY IMPLIES THAT DATA(1)=DATA(N/2+1)=0. ISIGN = +1 OR -1,
C THE DIRECTION OF THE TRANSFORM. THE FIRST N/2+1 POINTS OF THE
C TRANSFORM, REAL AND SAME SYMMETRY AS THE INPUT, WILL BE RETURNED
C IN DATA, REPLACING THE INPUT. RUNNING TIME AND ARRAY
C STORAGE ARE EQUIVALENT TO THE TRANSFORM OF N/4 COMPLEX, ASYM-
C METRIC DATA. NOTE THAT THE COMPLEX TRANSFORM OF EVEN OR ODD
C SYMMETRIC DATA IS EQUIVALENT TO THE COSINE OR SINE TRANSFORM,
C RESPECTIVELY, OF ASYMMETRIC DATA.
C
SUBROUTINE FOUR2 (DATA,N,NDIM,ISIGN,IFORM)
DOUBLE PRECISION DATA(*)
INTEGER N(*) !Read-only
INTEGER NDIM !Read-only
INTEGER ISIGN !Read-only
INTEGER IFORM !Read-only
C COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN.
C MULTI-DIMENSIONAL TRANSFORM, EACH DIMENSION A POWER OF TWO,
C COMPLEX OR REAL DATA.
C TRANSFORM(K1,K2,...) = SUM(DATA(J1,J2,...)*EXP(ISIGN*2*PI*SQRT(-1)
C *((J1-1)*(K1-1)/N(1)+(J2-1)*(K2-1)/N(2)+...))), SUMMED FOR ALL
C J1 AND K1 FROM 1 TO N(1), J2 AND K2 FROM 1 TO N(2),
C ETC. FOR ALL NDIM SUBSCRIPTS. NDIM MUST BE POSITIVE AND
C EACH N(IDIM) MUST BE A POWER OF TWO. ISIGN IS +1 OR -1.
C LET NTOT = N(1)*N(2)*...*N(NDIM). THEN A -1 TRANSFORM
C FOLLOWED BY A +1 ONE (OR VICE VERSA) RETURNS NTOT
C TIMES THE ORIGINAL DATA. IFORM = 1, 0 OR -1, AS DATA IS
C COMPLEX, REAL OR THE FIRST HALF OF A COMPLEX ARRAY. TRANSFORM
C VALUES ARE RETURNED TO ARRAY DATA. THEY ARE COMPLEX, REAL OR
C THE FIRST HALF OF A COMPLEX ARRAY, AS IFORM = 1, -1 OR 0.
C THE TRANSFORM OF A REAL ARRAY (IFORM = 0) DIMENSIONED N(1) BY N(2)
C BY ... WILL BE RETURNED IN THE SAME ARRAY, NOW CONSIDERED TO
C BE COMPLEX OF DIMENSIONS N(1)/2+1 BY N(2) BY .... NOTE THAT IF
C IFORM = 0 OR -1, N(1) MUST BE EVEN, AND ENOUGH ROOM MUST BE
C RESERVED. THE MISSING VALUES MAY BE OBTAINED BY COMPLEX CONJUGA-
C TION. THE REVERSE TRANSFORMATION, OF A HALF COMPLEX ARRAY DIMEN-
C SIONED N(1)/2+1 BY N(2) BY ..., IS ACCOMPLISHED BY SETTING IFORM
C TO -1. IN THE N ARRAY, N(1) MUST BE THE TRUE N(1), NOT N(1)/2+1.
C THE TRANSFORM WILL BE REAL AND RETURNED TO THE INPUT ARRAY.
C RUNNING TIME IS PROPORTIONAL TO NTOT*LOG2(NTOT), RATHER THAN
C THE NAIVE NTOT**2. FURTHERMORE, LESS ERROR IS BUILT UP.
C WRITTEN BY NORMAN BRENNER, MIT, JUNE 1968. SEE-- IEEE AUDIO
C TRANSACTIONS, JUNE 1967 AND JUNE 1969, SPECIAL FFT ISSUES.
C
SUBROUTINE BITRV (DATA,NPREV,N,NREM)
DOUBLE PRECISION DATA(*)
INTEGER NPREV !Read-only
INTEGER N !Read-only
INTEGER NREM !Read-only
C SHUFFLE THE DATA BY BIT REVERSAL.
C DIMENSION DATA(NPREV,N,NREM)
C COMPLEX DATA
C EXCHANGE DATA(J1,J4REV,J5) WITH DATA(J1,J4,J5) FOR ALL J1 FROM 1
C TO NPREV, ALL J4 FROM 1 TO N (WHICH MUST BE A POWER OF TWO), AND
C ALL J5 FROM 1 TO NREM. J4REV-1 IS THE BIT REVERSAL OF J4-1. E.G.
C SUPPOSE N = 32. THEN FOR J4-1 = 10011, J4REV-1 = 11001, ETC.
C
SUBROUTINE COOL2 (DATA,NPREV,N,NREM,ISIGN)
DOUBLE PRECISION DATA(*)
INTEGER NPREV !Read-only
INTEGER N !Read-only
INTEGER NREM !Read-only
INTEGER ISIGN !Read-only
C DISCRETE FOURIER TRANSFORM OF LENGTH N. IN-PLACE COOLEY-TUKEY
C ALGORITHM, BIT-REVERSED TO NORMAL ORDER, SANDE-TUKEY PHASE SHIFTS.
C DIMENSION DATA(NPREV,N,NREM)
C COMPLEX DATA
C DATA(J1,K4,J5) = SUM(DATA(J1,J4,J5)*EXP(ISIGN*2*PI*I*(J4-1)*
C (K4-1)/N)), SUMMED OVER J4 = 1 TO N FOR ALL J1 FROM 1 TO NPREV,
C K4 FROM 1 TO N AND J5 FROM 1 TO NREM. N MUST BE A POWER OF TWO.
C METHOD--LET IPREV TAKE THE VALUES 1, 2 OR 4, 4 OR 8, ..., N/16,
C N/4, N. THE CHOICE BETWEEN 2 OR 4, ETC., DEPENDS ON WHETHER N IS
C A POWER OF FOUR. DEFINE IFACT = 2 OR 4, THE NEXT FACTOR THAT
C IPREV MUST TAKE, AND IREM = N/(IFACT*IPREV). THEN--
C DIMENSION DATA(NPREV,IPREV,IFACT,IREM,NREM)
C COMPLEX DATA
C DATA(J1,J2,K3,J4,J5) = SUM(DATA(J1,J2,J3,J4,J5)*EXP(ISIGN*2*PI*I*
C (K3-1)*((J3-1)/IFACT+(J2-1)/(IFACT*IPREV)))), SUMMED OVER J3 = 1
C TO IFACT FOR ALL J1 FROM 1 TO NPREV, J2 FROM 1 TO IPREV, K3 FROM
C 1 TO IFACT, J4 FROM 1 TO IREM AND J5 FROM 1 TO NREM. THIS IS
C A PHASE-SHIFTED DISCRETE FOURIER TRANSFORM OF LENGTH IFACT.
C FACTORING N BY FOURS SAVES ABOUT TWENTY FIVE PERCENT OVER FACTOR-
C ING BY TWOS. DATA MUST BE BIT-REVERSED INITIALLY.
C IT IS NOT NECESSARY TO REWRITE THIS SUBROUTINE INTO COMPLEX
C NOTATION SO LONG AS THE FORTRAN COMPILER USED STORES REAL AND
C IMAGINARY PARTS IN ADJACENT STORAGE LOCATIONS. IT MUST ALSO
C STORE ARRAYS WITH THE FIRST SUBSCRIPT INCREASING FASTEST.
C
SUBROUTINE FIXRL (DATA,N,NREM,ISIGN,IFORM)
DOUBLE PRECISION DATA(*)
INTEGER N !Read-only
INTEGER NREM !Read-only
INTEGER ISIGN !Read-only
INTEGER IFORM !Read-only
C FOR IFORM = 0, CONVERT THE TRANSFORM OF A DOUBLED-UP REAL ARRAY,
C CONSIDERED COMPLEX, INTO ITS TRUE TRANSFORM. SUPPLY ONLY THE
C FIRST HALF OF THE COMPLEX TRANSFORM, AS THE SECOND HALF HAS
C CONJUGATE SYMMETRY. FOR IFORM = -1, CONVERT THE FIRST HALF
C OF THE TRUE TRANSFORM INTO THE TRANSFORM OF A DOUBLED-UP REAL
C ARRAY. N MUST BE EVEN.
C USING COMPLEX NOTATION AND SUBSCRIPTS STARTING AT ZERO, THE
C TRANSFORMATION IS--
C DIMENSION DATA(N/2,NREM)
C ZSTP = EXP(ISIGN*2*PI*I/N)
C DO 10 I2=0,NREM-1
C DATA(0,I2) = CONJ(DATA(0,I2))*(1+I)
C DO 10 I1=1,N/4
C Z = (1+(2*IFORM+1)*I*ZSTP**I1)/2
C I1CNJ = N/2-I1
C DIF = DATA(I1,I2)-CONJ(DATA(I1CNJ,I2))
C TEMP = Z*DIF
C DATA(I1,I2) = (DATA(I1,I2)-TEMP)*(1-IFORM)
C 10 DATA(I1CNJ,I2) = (DATA(I1CNJ,I2)+CONJ(TEMP))*(1-IFORM)
C IF I1=I1CNJ, THE CALCULATION FOR THAT VALUE COLLAPSES INTO
C A SIMPLE CONJUGATION OF DATA(I1,I2).
C
SUBROUTINE FIT1 (NIT,OMEGA,ACF,MS,NTERMS,W,A00,SIGA00,NSC,K,QSSW,YFIT,A,SIGA,
COVMAT,CORMAT,CHISQR,QGOOD)
INTEGER NIT
DOUBLE PRECISION OMEGA(*)
DOUBLE PRECISION ACF(*) !Read-only
INTEGER MS !Read-only
INTEGER NTERMS !Read-only
DOUBLE PRECISION W(*) !Read-only
DOUBLE PRECISION A00(NFPMAX) !Read-only
DOUBLE PRECISION SIGA00(NFPMAX) !Read-only
INTEGER NSC(*) !Read-only
INTEGER K !Read-only
INTEGER QSSW(*) !Read-only
DOUBLE PRECISION YFIT(*)
DOUBLE PRECISION A(NFPMAX)
DOUBLE PRECISION SIGA(NFPMAX)
DOUBLE PRECISION COVMAT(NFPMAX*NFPMAX)
DOUBLE PRECISION CORMAT(NFPMAX*NFPMAX)
DOUBLE PRECISION CHISQR
LOGICAL QGOOD
C
C JMH - 9/82 HARRIS FORTRAN 77
C
C FIT1 COMPUTES A NON-LINEAR-LEAST-SQUARES FIT TO A MEASURED ACF.
C
C NIT - MAXIMUM NUMBER OF ITERATIONS (INPUT)
C ACTUAL NUMBER OF ITERATIONS (OUTPUT)
C OMEGA - FREQUENCIES OF MEASURED ACF
C ACF - MEASURED ACF
C MS - NUMBER OF PAIRS OF DATA POINTS (ACF)
C NTERMS - NUMBER OF SEARCH PARAMETERS
C W - ACF WEIGHTS
C A00 - INITIAL VALUES OF FIT PARAMETERS
C SIGA00 - INITIAL VALUES OF FIT PARAMETER STANDARD DEVIATIONS
C NSC - NSC(I)=1 -> SEARCH PARAMETER A(I)
C K - REQUIRED LEVEL OF ACCURACY OF FIT CONVERGENCE (1-4)
C QSSW - INSCAL CONTROL SWITCH ARRAY
C
C YFIT - ACF FIT
C A - OUTPUT FIT PARAMETERS
C SIGA - OUTPUT STANDARD DEVIATIONS
C COVMAT - COVARIANCE MATRIX
C CORMAT - CORRELATION MATRIX
C CHISQR - FIT CHISQR
C QGOOD - 1 - FIT WAS SUCCESSFUL
C 2 - FIT WAS UNSUCCESSFUL
C
SUBROUTINE FITPH (ACFR,ACFI,WR,WI,NLAG,DTAU,V,DV)
DOUBLE PRECISION ACFR(*) !Read-only
DOUBLE PRECISION ACFI(*) !Read-only
DOUBLE PRECISION WR(*) !Read-only
DOUBLE PRECISION WI(*) !Read-only
INTEGER NLAG !Read-only
DOUBLE PRECISION DTAU !Read-only
DOUBLE PRECISION V
DOUBLE PRECISION DV
C
C JMH - 3/99
C
SUBROUTINE FITSTATS (SNLIM,CHISQR,CHISQT)
DOUBLE PRECISION SNLIM !Read-only
DOUBLE PRECISION CHISQR(*) !Read-only
DOUBLE PRECISION CHISQT(*) !Read-only
C
C JMH - 9/97
C
SUBROUTINE FIXIFAM (ASET)
DOUBLE PRECISION ASET(NFPMAX,6) !Read-only
C
C FIXIFAM APPLIES CORRECTIONS TO THE IFAM STRUCTURE. IDEALLY,
C THIS SHOULD NOT BE NECESSARY. FOR EXAMPLE, IF IFAM IS PRODUCED
C BY CMST2IFAM, IN MOST CASES CMST2IFAM SHOULD BE MODIFIED TO
C CREATE A GOOD IFAM. PRACTICALLY SPEAKING, THIS IS NOT ALWAYS
C POSSIBLE. FIXIFAM IS IN SOME SENSE A KLUDGE ROUTINE, AND
C MAY HAVE TO BE MODIFIED WHEN DATA FROM A NEW SOURCE IS
C ANALYZED. A PRIMARY VIRTUE IS THAT IT COLLECTS SEVERAL KLUDGES
C TOGETHER WHERE THEY ARE MORE EASILY MONITORED THAN IF THEY WERE
C SPREAD THROUGHOUT INSCAL.
C
C JMH - 8/97
C JMH - 12/98 - GENERALIZED FOR INSCAL 8.0
C
SUBROUTINE FIXSPEC (ACF,NLAG,DTAU,IND,ASYM,SPECFAC,NSPECT,OMEGA,SPECT,AACF)
DOUBLE PRECISION ACF(2,*) !Read-only
INTEGER NLAG !Read-only
DOUBLE PRECISION DTAU !Read-only
INTEGER IND !Read-only
DOUBLE PRECISION ASYM !Read-only
DOUBLE PRECISION SPECFAC(*) !Read-only
INTEGER NSPECT
DOUBLE PRECISION OMEGA(*)
DOUBLE PRECISION SPECT(2,*)
DOUBLE PRECISION AACF(2,*)
C
C JMH - 8/97
C
C SPECTCAL COMPUTES THE SPECTRUM CORRESPONDING TO A COMPLEX,
C HERMITIAN AUTOCORRELATION FUNCTION SPECIFIED BY ITS POSITIVE
C LAGS.
C
C INPUT:
C ACF - INPUT AUTOCORRELATION FUNCTION. ACF(1,I) CONTAINS
C THE REAL PART OF THE ACF AND ACF(2,I) CONTAINS
C THE IMAGINARY PART OF THE ACF.
C NLAG - NUMBER OF LAGS IN ACF. IF NLAG IS NOT A
C POWER OF 2 THE ACF WILL BE ZERO PADDED TO
C THE NEXT HIGHER POWER OF 2 (=NLAGP2).
C DTAU - LAG SPACING (SEC)
C IND=0 - NO CORRECTION
C 1 - APPLY ASSYMETRY FACTOR
C 2 - APPLY SPECTRAL CORRECTION FACTOR
C ASYM - SPECTRUM ASSYMETRY FACTOR.
C
C OUTPUT:
C NSPECT - NUMBER OF POINTS IN THE OUTPUT SPECTRUM. IF
C NLAG IS A POWER OF 2, NSPECT=2*NLAG. OTHERWISE,
C NSPECT = 2*NLAGP2.
C OMEGA - ANGULAR FREQUENCY (OMEGA=2*PI*FREQ).
C SPECT - OUTPUT SPECTRUM. THE SECOND DIMENSION OF
C SPECT MUST BE AT LEAST EQUAL TO NSPECT.
C AACF - OUTPUT AUTOCORRELATION FUNCTION. ACF(1,I) CONTAINS
C THE REAL PART OF THE ACF AND ACF(2,I) CONTAINS
C THE IMAGINARY PART OF THE ACF.
C
C DO 100 I=1,NLAG
C WRITE(30,"(1X,I4,3E13.5)") I, (I-1)*DTAU, ACF(1,I), ACF(2,I)
C 100 CONTINUE
C
SUBROUTINE FOUR1 (DATA,NN,ISIGN)
DOUBLE PRECISION DATA(2*NN)
INTEGER NN !Read-only
INTEGER ISIGN !Read-only
SUBROUTINE FULLSPEC (ACF,NLAG,DTAU,ASYM,NSPECT,OMEGA,SPECT,AACF)
DOUBLE PRECISION ACF(2,*) !Read-only
INTEGER NLAG !Read-only
DOUBLE PRECISION DTAU !Read-only
DOUBLE PRECISION ASYM !Read-only
INTEGER NSPECT
DOUBLE PRECISION OMEGA(*)
DOUBLE PRECISION SPECT(2,*)
DOUBLE PRECISION AACF(2,*)
C
C JMH - 8/97
C
C SPECTCAL COMPUTES THE SPECTRUM CORRESPONDING TO A COMPLEX,
C HERMITIAN AUTOCORRELATION FUNCTION SPECIFIED BY ITS POSITIVE
C LAGS.
C
C INPUT:
C ACF - INPUT AUTOCORRELATION FUNCTION. ACF(1,I) CONTAINS
C THE REAL PART OF THE ACF AND ACF(2,I) CONTAINS
C THE IMAGINARY PART OF THE ACF.
C NLAG - NUMBER OF LAGS IN ACF. IF NLAG IS NOT A
C POWER OF 2 THE ACF WILL BE ZERO PADDED TO
C THE NEXT HIGHER POWER OF 2 (=NLAGP2).
C DTAU - LAG SPACING (SEC)
C ASYM - SPECTRUM ASSYMETRY FACTOR.
C
C OUTPUT:
C NSPECT - NUMBER OF POINTS IN THE OUTPUT SPECTRUM. IF
C NLAG IS A POWER OF 2, NSPECT=2*NLAG. OTHERWISE,
C NSPECT = 2*NLAGP2.
C OMEGA - ANGULAR FREQUENCY (OMEGA=2*PI*FREQ).
C SPECT - OUTPUT SPECTRUM. THE SECOND DIMENSION OF
C SPECT MUST BE AT LEAST EQUAL TO NSPECT.
C AACF - OUTPUT AUTOCORRELATION FUNCTION. ACF(1,I) CONTAINS
C THE REAL PART OF THE ACF AND ACF(2,I) CONTAINS
C THE IMAGINARY PART OF THE ACF.
C
C DO 100 I=1,NLAG
C WRITE(30,"(1X,I4,3E13.5)") I, (I-1)*DTAU, ACF(1,I), ACF(2,I)
C 100 CONTINUE
C
SUBROUTINE FUNCTLWE (NPTS,X,NTYPE,NTERMS,A,NS,FIT)
INTEGER NPTS !Read-only
DOUBLE PRECISION X(*) !Read-only
INTEGER NTYPE(*) !Read-only
INTEGER NTERMS !Read-only
DOUBLE PRECISION A(*) !Read-only
INTEGER NS(*) !Read-only
DOUBLE PRECISION FIT(*)
C
C JMH - 11/93
C
C FUNCTLWE CALCULATES THE CURFIT AUTOCORRELATION OBJECT FUNCTION.
C
SUBROUTINE FUNCTN (MS,OMEGA,NTYPE,NTERMS,A,NS,ACF)
INTEGER MS !Read-only
DOUBLE PRECISION OMEGA(*)
INTEGER NTYPE(*) !Read-only
INTEGER NTERMS !Read-only
DOUBLE PRECISION A(*) !Read-only
INTEGER NS(*) !Read-only
DOUBLE PRECISION ACF(*)
C
C JMH - 9/81
C JMH - 4/92 MIDAS Version
C
C FUNCTN CALCULATES THE CURFIT AUTOCORRELATION OBJECT FUNCTION.
C
C KIND - INDICATES SEARCH PARAMETERS, SPECTRUM CALCULATION
C FUNCTION
C 1 - M+
C 2 - M+, O+
C 3 - O+
C 4 - O+, H+
C 5 - M+ COLLISIONS INCLUDED
C JIND - INDICATES WHETHER VELOCITY IS TO BE SEARCHED
C 0 - V(O+)=V(H+)=0
C 1 - V(O+)
C
SUBROUTINE FUNPH (NTAU,TAU,NTYPE,NTERMS,A,NS,PHASE)
INTEGER NTAU !Read-only
DOUBLE PRECISION TAU(*) !Read-only
INTEGER NTYPE(*) !Read-only
INTEGER NTERMS !Read-only
DOUBLE PRECISION A(*) !Read-only
INTEGER NS(*) !Read-only
DOUBLE PRECISION PHASE(*)
C
C JMH - 3/99
C
C FUNPH CALCULATES THE CURFIT AUTOCORRELATION PHASE OBJECT FUNCTION.
C
SUBROUTINE GEOPRM (IYDI,SECI,GLATI,GLONGI)
INTEGER IYDI !Read-only
DOUBLE PRECISION SECI !Read-only
DOUBLE PRECISION GLATI !Read-only
DOUBLE PRECISION GLONGI !Read-only
C
SUBROUTINE GTCMTM (JDAYNC,STRSEC,ENDSEC)
INTEGER JDAYNC !Unknown (program incomplete)
INTEGER STRSEC
INTEGER ENDSEC
C
C GTCMTM COMPUTES THE START AND END TIME OF DATA TO BE ANALYZED
C BY INSCAL. THESE ARE CALCULATED FROM THE START DATE, START TIME
C AND END TIME SPECIFIED ON THE INSCAL COMMAND LINE OR FILE.
C
C JDAYNC - START JULIAN DAY NUMBER
C STRSEC - START TIME IN SECONDS RELATIVE TO MIDNIGHT ON JDAYNC
C ENDSEC - END TIME IN SECONDS RELATIVE TO MIDNIGHT ON JDAYNC
C
SUBROUTINE GTIFTM (JDAYNC,SEC)
INTEGER JDAYNC !Read-only
INTEGER SEC
C
SUBROUTINE GTMDAT (MDFILE,ISTAT)
CHARACTER MDFILE !Read-only
INTEGER ISTAT
C
SUBROUTINE GETEXT (TXTFIL,CODE,INDEX,TEXT)
CHARACTER * 80 TXTFIL !Read-only
CHARACTER * 16 CODE(*)
INTEGER INDEX(*)
CHARACTER * 128 TEXT(*)
C
SUBROUTINE PRNTEXT (LFN,IC,CODE,INDEX,TEXT)
INTEGER LFN !Read-only
CHARACTER IC !Unknown (program incomplete)
CHARACTER * 16 CODE(*) !Unknown (program incomplete)
INTEGER INDEX(*) !Read-only
CHARACTER * 128 TEXT(*) !Unknown (program incomplete)
C
SUBROUTINE IF2OFI (NZR,Z1,Z2)
INTEGER NZR !Read-only
DOUBLE PRECISION Z1(*) !Read-only
DOUBLE PRECISION Z2(*) !Read-only
C
C IF2OFI INITIALIZES OFAM FROM IFAM. ALL PARAMETERS WHICH CAN
C BE PASSED DIRECTLY FROM IFAM TO OFAM ARE HANDLED HERE.
C PARAMETERS WHICH DEPEND ON ANALYSIS OF THE DATA IN IFAM
C ARE HANDLED ELSEWHERE.
C
SUBROUTINE IFAMINF (INFILE)
CHARACTER INFILE !Read-only
C
C JMH - 5/13/91
C
C IFAMINF READS AN ASCII FORMAT IFAM FILE INTO
C LABELLED COMMON IFAM.
C
SUBROUTINE IFAMOUTF (OTFILE)
CHARACTER OTFILE !Read-only
C
C JMH - 5/13/91
C
C IFAMOUTF WRITES LABELLED COMMON IFAM TO AN ASCII FORMAT
C IFAM FILE.
C
SUBROUTINE IFAMOUTNF (LFN)
INTEGER LFN !Read-only
C
C JMH - 5/13/91
C
C IFAMOUTF WRITES LABELLED COMMON IFAM TO AN ASCII FORMAT
C IFAM FILE.
C
SUBROUTINE INICOM ()
C
C JMH - 10/83
C JMH - 4/92 RENAMED FROM BLKDAT TO INICOM
C
C INICOM INITIALIZES THE VARIABLES IN /DATA/ /STAT/, /QSSWEQ/ and
C /PARCOM/.
C
SUBROUTINE KNOTSION (Z1,DZ1,Z2,DZ2,NDIM,N,Z)
DOUBLE PRECISION Z1 !Read-only
DOUBLE PRECISION DZ1 !Read-only
DOUBLE PRECISION Z2 !Read-only
DOUBLE PRECISION DZ2 !Read-only
INTEGER NDIM !Read-only
INTEGER N
DOUBLE PRECISION Z(NDIM)
C
C JMH - 01/2004
C
C KNOTSION CALCULATES A SET OF SPLINE KNOTS SUITABLE FOR
C FITTING IONOSPHERIC PARAMETER PROFILES.
C
C INPUT PARAMETERS
C A,B - EXPONENTIAL KNOT-SPACING PARAMETERS
C DZ = 10.0D0**(A+B*Z)
C A=0.8, B=0.002 -> DZ=10 KM AT 100 KM,
C DZ=100 KM AT 600 KM
C Z1 - LOWER ALTITUDE BOUND (FIRST KNOT)
C Z2 - UPPER ALTITUDE BOUND (LAST KNOT)
C NDIM - MAXIMUM NUMBER OF KNOTS
C
C OUTPUT PARAMETERS
C N - NUMBER OF KNOTS
C Z - KNOTS
C
SUBROUTINE FTL1L2 (MDIM,NDIM,M,N,X,Y,DY,BASIS,NB,CB,IFTYPE,RESLIM,YFIT,A,
COVMAT,RNORM,IND,IWORK,WORK)
INTEGER MDIM !Read-only
INTEGER NDIM !Read-only
INTEGER M !Read-only
INTEGER N !Read-only
DOUBLE PRECISION X(MDIM) !Read-only
DOUBLE PRECISION Y(MDIM) !Read-only
DOUBLE PRECISION DY(MDIM) !Read-only
SUBROUTINE BASIS
Argument 1: INTEGER
Argument 2: INTEGER
Argument 3: DOUBLE PRECISION Array Element
Argument 4: DOUBLE PRECISION Array Element
INTEGER NB !Read-only
DOUBLE PRECISION CB(MDIM,*) !Read-only
INTEGER IFTYPE !Read-only
DOUBLE PRECISION RESLIM !Read-only
DOUBLE PRECISION YFIT(MDIM)
DOUBLE PRECISION A(NDIM)
DOUBLE PRECISION COVMAT(NDIM,NDIM)
DOUBLE PRECISION RNORM
INTEGER IND
INTEGER IWORK(*)
DOUBLE PRECISION WORK(*)
C
C JMH - 9/88
C
C FTL1L2 COMPUTES A WEIGHTED LINEAR L1 (LEAST ABSOLUTE VALUES) OR
C L2 (LEAST SQUARES) FIT TO AN ARRAY OF DATA. A ROBUST L2 FIT
C MAY BE COMPUTED BY USING THE L1 FIT TO REJECT DATA OUTLIERS.
C
C INPUT PARAMETERS
C MDIM - MAXIMUM NUMBER OF DATA POINTS AS DIMENSIONED
C IN THE CALLING PROGRAM
C NDIM - MAXIMUM NUMBER OF FIT COEFFICIENTS AS DIMENSIONED
C IN THE CALLING PROGRAM
C M - NUMBER OF DATA POINTS
C N - NUMBER OF FIT COEFFICIENTS TO BE COMPUTED
C X - ARRAY OF INDEPENDENT VARIABLE VALUES (USED IN BASIS
C FUNCTION CALCULATIONS
C Y - ARRAY OF DEPENDENT VARIABLE VALUES (THE DATA TO BE FIT
C DY - ARRAY OF DATA ERRORS, USED TO WEIGHT THE FITS
C BASIS - BASIS FUNCTION ROUTINE BASIS(N,I,X,F)
C N - NUMBER OF BASIS FUNCTIONS TO CALCULATE (INPUT)
C I - INDEX INTO X,F ARRAYS (INPUT)
C X - X VALUE AT WHICH TO EVALUATE BASIS FUNCTION (IN
C F - ARRAY OF N BASIS FUNCTION VALUES (OUTPUT)
C IN SIMPLE CASES, SUCH AS A POLYNOMIAL IN X, BASIS MAY
C ABLE TO EVALUATE F GIVEN ONLY X. IN MORE COMPLEX CASE
C ADDITIONAL INFORMATION MAY BE REQUIRED. IN THESE CASES
C A COMMON AREA MAY BE USED TO PASS THE REQUIRED INFORMA
C FROM THE ROUTINE WHICH CALLS FTL1L2 DIRECTLY TO BASIS.
C THESE CASES I MAY BE NEEDED TO CORRECTLY ACCESS TARRAY
C THE COMMON BLOCK.
C NB - NUMBER OF BASIS VECTORS
C CB - COEFFICIENTS OF BASIS VECTORS IN BASIS FUNCTION.
C I.E., SPLINE FUNCTION BASIS FUNCTIONS CAN RETURN UP TO
C FOUR BASIS VECTORS, THE 0TH THROUGH THIRD DERIVATIVES.
C THE DERIVATIVES ARE OFTEN USED TO CONSTRAIN THE SPLINE.
C MORE GENERALLY, CB ALLOW DATA WHICH ARE A LINEAR
C COMBINATION OF DERIVATIVES TO BE FIT, AS MIGHT BE
C REQUIRED TO CONSTRAIN THE FIT TO OBEY A DIFFERENTIAL
C EQUATION.
C IFTYPE - 1 -> DO L1 FIT
C 2 -> DO L2 FIT
C 3 -> DO L1 FIT, REJECT DATA FARTHER THAN RESLIM FROM
C THE MEDIAN L1 RESIDUAL, DO L2 FIT TO REMAINING DA
C RESLIM - RESIDUAL LIMIT FOR ACCEPTING DATA FOR ROBUST L2 FIT
C
C OUTPUT PARAMETERS
C YFIT - ARRAY OF FIT VALUES (THE BEST FIT TO Y)
C A - ARRAY OF FIT COEFFICIENTS (THE SOLUTION VECTOR)
C COVMAT - COVARIANCE MATRIX OF L2 FIT
C RNORM - L1 NORM (IFTYPE=1) OR L2 NORM (IFTYPE=2,3) OF THE
C RESIDUAL VECTOR
C IND - 10*IND1 + IND2
C IND1= 0 -> L1 SOLUTION VECTOR IS UNIQUE
C =+1 -> L1 SOLUTION VECTOR PROBABLY NOT UNIQUE
C =-1 -> L1 FIT TERMINATED PREMATURELY
C IND2= 0 -> L2 FIT SUCCESSFUL
C
C WORK ARRAYS
C IWORK - DIMENSIONED AT LEAST 2MDIM*NDIM
C WORK - DIMENSIONED AT LEAST NDIM*MDIM+4*NDIM+5*MDIM+NDIM**2+
C MDIM*(MDIM+3)/2+4*MDIM
C
SUBROUTINE SHELL1 (A,N)
DOUBLE PRECISION A(*)
INTEGER N !Read-only
C
C SHELL1 IS AN ADAPTATION OF THE SORTING ROUTINE SHELLSORT
C (ACM ALGORITHM 201) BY J. BOOTHROYD. THE CONTENTS OF THE ARRAY
C A, COMPRISING N ELEMENTS IN N SUCCESSIVE WORDS, WILL BE SORTED
C INTO NON-DESCENDING ORDER.
C
SUBROUTINE PSLV (ID,K,MM,MMM,P,B,X)
INTEGER ID !Read-only
INTEGER K !Read-only
INTEGER MM !Read-only
INTEGER MMM !Read-only
DOUBLE PRECISION P(MMM) !Read-only
DOUBLE PRECISION B(MM) !Read-only
DOUBLE PRECISION X(MM)
C this subroutine solves the square non-singular system
C of linear equations
C p*x=b,
C or the square non-singular system of linear equations
C p(transpose)*x=b,
C where p is an upper triangular matrix, b is the right hand
C side vector and x is the solution vector.
C the input data to the subroutine.
C id an integer indicator specified by the user.
C if id=1 the equation p*x=b is solved.
C if id= any integer other than 1, the equation
C p(transpose)*x=b is solved.
C k an integer = the number of equations of the given
C system.
C mm an integer greater than or equal to k.
C mmm an integer = (mm*(mm+3))/2
C p an mmm-vector. its first (k+1) elements contain the
C first k elements of row 1 of the upper triangular
C matrix plus an extra location to the right. its
C next k elements contain the (k-1) elements of row 2
C of the upper triangular matrix plus an extra
C location to the right, ..., etc.
C b an mm-vector. its first k elements contain the
C r.h.s. vector of the given system.
C the output of the subroutine.
C x an mm-vector. on exit, its first k elements contain
C the solution to the given system.
C
SUBROUTINE L1 (M,N,MM,MMM,NN,CT,F,PREC,EPS,IC,IR,IB,UF,BP,XP,T,ALFA,P,GINV,
VB,IRANK,ITER,R,A,Z,IND)
INTEGER M !Read-only
INTEGER N !Read-only
INTEGER MM !Read-only
INTEGER MMM !Read-only
INTEGER NN !Read-only
DOUBLE PRECISION CT(MM,NN)
DOUBLE PRECISION F(NN) !Read-only
DOUBLE PRECISION PREC !Read-only
DOUBLE PRECISION EPS !Read-only
INTEGER IC(NN)
INTEGER IR(MM)
INTEGER IB(NN)
DOUBLE PRECISION UF(MM)
DOUBLE PRECISION BP(MM)
DOUBLE PRECISION XP(MM)
DOUBLE PRECISION T(NN)
DOUBLE PRECISION ALFA(NN)
DOUBLE PRECISION P(MMM)
DOUBLE PRECISION GINV(MM,MM)
DOUBLE PRECISION VB(MM)
INTEGER IRANK
INTEGER ITER
DOUBLE PRECISION R(NN)
DOUBLE PRECISION A(MM)
DOUBLE PRECISION Z
INTEGER IND
C this subroutine solves an overdetermined system of
C linear equations in the l1 norm by using a dual simplex
C algorithm to the linear programming formulation of the
C given problem. in this algorithm certain intermediate
C simplex iterations are skipped.
C for purpose of numerical stability, this subroutine
C uses a triangular decomposition to the basis matrix.
C the system of linear equations has the form
C c*a=f,
C where c is a given real n by m matrix of rank k.le.m.le.n
C and f is a given real n-vector.
C the problem to be solved is to calculate the elements
C of the m-vector a* which gives the minimum norm z.
C z= abs(r(1)) + abs(r(2)) + ... + abs(r(n)) ,
C where r(i) is the ith residual and is given by
C r(i)=c(i,1)*a(1)+c(i,2)*a(2)+ ... +c(i,m)*a(m)-f(i).
C subroutine l1 is completely self contained (consists
C of two subroutines l1 and pslv).
C the input data to the subroutine.
C m the number of columns of matrix c.
C n the number of rows of matrix c.
C mm an integer greater than or equal to m.
C mmm an integer = (mm*(mm+3))/2 .
C nn an integer greater than or equal to n.
C ct a matrix of dimensions mm by nn. on entry, its first
C m rows and first n columns contain the transpose of
C matrix c in the given system c*a=f.
C f an nn-vector. on entry, its first n elements contain
C the r.h.s. of the given system c*a=f.
C prec the round-off level of the computer. for the ibm
C 360/67 computer, prec is about 1.e-6 and 1.e-16 for
C single and double precision respectively.
C eps a specified tolerance such that a calculated number
C x is considered zero if abs(x).LT.eps. for the ibm 360/67
C computer, eps is usually taken 1.e-4 and 1.e-11
C respectively.
C the results of the problem.
C irank the calculated rank of matrix c.
C iter the number of iterations which the solution needed.
C a an mm-vector. its first m elements are the solution
C vector a*.
C r an nn-vector. its first n elements contain the
C the residual vector r=c*a-f.
C z the minimum l1 norm of the residual vector r.
C ind a return indicator. ind=0 indicates that the
C solution vector a* is unique. ind=1, indicates that
C a* is most probably not unique. ind=-1 indicates
C premature termination of the calculation due to
C very ill-conditioning of matrix c.
C the meaning of the other parameters.
C ginv an mm-square matrix. its first irank columns and
C first irank rows contain the inverse of the initial
C basis matrix and its update.
C vb an mm-vector. its first irank elements contain the
C initial basic solution and its update.
C t an nn-vector. its first n elements contain the
C elements of the row in the simplex tableau, that
C corresponds to the column which leaves the basis.
C alfa an nn-vector. its first n elements contain the
C ratios: alfa(j) = r(j)/t(j).
C ic an nn-vector. its first n elements contain the
C column indices of matrix ct.
C ir an mm vector. its first irank elements contain the
C row indices of the linearly independent rows of ct.
C ib a sign nn-vector. its first n elements have the
C values +1 or -1. ib(j)=+1 indicates that column j
C in the tableau is at its lower bound 0. ib(j)=-1
C indicates that column j is at its upper bound 2.
C p an mmm-vector. its first ((irank*(irank+3))/2)-1
C elements contain the (irank*(irank+1))/2 elements of
C the upper triangular matrix + extra (irank-1) working
C locations. see the comments in subroutine pslv.
C bp an mm-vector. its first irank elements are the
C r.h.s. of the triangular equations as p*xp=bp.
C xp an mm-vector . its first irank elements are the
C solution of the triangular equations as p*xp=bp.
C uf an mm-working vector.
SUBROUTINE INCLDC (N,KLEFT,KRIGHT,KBAND,KBIAS,WVAL,X,YVAL,RBAR,SS)
INTEGER N !Read-only
INTEGER KLEFT !Read-only
INTEGER KRIGHT !Read-only
INTEGER KBAND !Read-only
INTEGER KBIAS !Read-only
DOUBLE PRECISION WVAL !Read-only
DOUBLE PRECISION X(*)
DOUBLE PRECISION YVAL !Read-only
DOUBLE PRECISION RBAR(*)
DOUBLE PRECISION SS
C
C JMH - 6/80 ANS FORTRAN 66
C 9/88 IMPROVED COMMENTS. CODE UNCHANGED.
C
C INCLDC USES GIVENS TRANSFORMS TO COMPUTE THE LEAST SQUARES SOLUTIO
C OF AN OVERDETERMINED LINEAR SYSTEM. INCLUDC MUST BE CALLED ONCE F
C EACH EQUATION IN THE SYSTEM, THAT IS DATA ARE INCLUDED IN THE FIT
C AT A TIME. AFTER ALL DATA HAVE BEEN INCLUDED, REGRSC MAY BE CALLE
C TO COMPUTE THE FIT COEFFICIENTS, AND COVC MAY BE CALLED TO COMPUTE
C THE FIT COVARIANCE MATRIX. THIS VERSION OF THE "INCLUDE" ALGORITHM
C PROVIDES FOR THE EFFICIENT HANDLING OF BANDED SYSTEMS. SPLBAS IS A
C B-SPLINE BASIS FUNCTION ROUTINE WHICH MAKES USE OF THIS FEATURE.
C
C INPUT PARAMETERS:
C N - NUMBER OF TERMS IN EACH EQUATION
C KLEFT - LEFTMOST NON-ZERO ELEMENT OF X (=1 FOR UNBANDED SYSTEM
C KRIGHT - RIGHTMOST NON-ZERO ELEMENT OF X (=N FOR UNBANDED SYSTE
C KBAND - BANDWIDTH OF THE SYSTEM (=N FOR UNBANDED SYSTEM)
C KBIAS - BIAS INTO X (USUALLY 0)
C WVAL - WEIGHT OF THE EQUATION (TYPICALLY 1/YVAL_ERROR**2)
C X - LEFT-HAND-SIDE OF EQUATION TO BE INCLUDED
C YVAL - RIGHT-HAND-SIDE OF EQUATION TO BE INCLUDED
C
C OUTPUT PARAMETERS:
C RBAR - AUGMENTED TRIANGULAR SYSTEM. MUST BE INITIALIZED TO
C ZERO BEFORE INCLDC IS CALLED FOR THE FIRST TIME.
C MUST BE DIMENSIONED AT LEAST (KBAND+2)*N-KBAND*(KBAND-
C SS - SUM-OF-SQUARES DEVIATION.
C
C IMPORTANT LOCAL VARIABLES:
C IT - IT<0 IF ROW HAS ZEROES. IT>0 IF ROW HAS NO ZEROES.
C ILZ - INDEX OF LAST ROW WITH ZEROES.
C NRZ0 - TOTAL NUMBER OF ELEMENTS IN ROWS WITH ZEROES.
C IB - MAGIC NUMBER USED TO COMPUTE INDEX (NT) OF LAST
C ELEMENT IN ROW.
C ND - INDEX OF FIRST ELEMENT IN ROW.
C NT - INDEX OF LAST ELEMENT IN ROW.
C I - ROW INDEX.
C
C EXAMPLE:
C
C RBAR I ND NT IT
C
C D R R R R 0 0 0 0 0 0 T 1 1 6 -5 N = 10
C D R R R R 0 0 0 0 0 T 2 7 12 -4 KBAND = 4
C D R R R R 0 0 0 0 T 3 13 18 -3
C D R R R R 0 0 0 T 4 19 24 -2
C D R R R R 0 0 T 5 25 30 -1
C D R R R R 0 T 6 31 36 0 NRZ0 = 36
C D R R R R T 7 37 42 1
C D R R R T 8 43 47 2
C D R R T 9 48 51 3
C D R T 10 52 54 4
C
SUBROUTINE SROTMG (SD1,SD2,SX1,SY1,SPARAM)
DOUBLE PRECISION SD1
DOUBLE PRECISION SD2
DOUBLE PRECISION SX1
DOUBLE PRECISION SY1 !Read-only
DOUBLE PRECISION SPARAM(5)
C
C CONSTRUCT THE MODIFIED GIVENS TRANSFORMATION MATRIX H WHICH ZEROS
C THE SECOND COMPONENT OF THE 2-VECTOR (SQRT(SD1)*SX1,SQRT(SD2)*
C SY1)**T).
C WITH SPARAM((1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
C
C SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0
C
C (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0)
C H=( ) ( ) ( ) ( )
C (SH21 SH22) (SH21 1.E0) (-1.E0 SH22) (0.E0 1.E0)
C LOCATIONS 2-5 OF SPARAM CONTAIN SH11, SH21, SH12 AND SH22
C RESPECTIVELY. (VALUES OF 1.E0, -1.E0, OR 0.E0 IMPLIED BY THE
C VALUE OF SPARAM(1) ARE NOT STORED IN SPARAM.)
C
C THE VALUES OF GAMSQ AND RGAMSQ SET IN THE DATA STATEMENT MAY BE
C INEXACT. THIS IS OK AS THEY ARE ONLY USED FOR TESTING THE SIZE
C OF SD1 AND SD2. ALL ACTUAL SCALING OF DATA IS DONE USING GAM.
C
SUBROUTINE SROTM (N,SX,INCX,SY,INCY,SPARAM)
INTEGER N !Read-only
DOUBLE PRECISION SX(*)
INTEGER INCX !Read-only
DOUBLE PRECISION SY(*)
INTEGER INCY !Read-only
DOUBLE PRECISION SPARAM(5) !Read-only
C
C APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
C
C (SX**T), WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF SX ARE IN
C (DX**T)
C
C SX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
C LX = (-INCX)*N, AND SIMILARLY FOR SY USING LX AND INCY.
C WITH SPARAM(1)=SFLAG, H HAS ONE OF THE FOLLOWING FORMS..
C
C SFLAG=-1.E0 SFLAG=0.E0 SFLAG=1.E0 SFLAG=-2.E0
C (SH11 SH12) (1.E0 SH12) (SH11 1.E0) (1.E0 0.E0)
C H=( ) ( ) ( ) ( )
C (SH21 SH22) (SH21 1.E0) (-1.E0 SH22) (0.E0 1.E0)
C SEE SROTMG FOR A DESCRIPTION OF DATA STORAGE IN SPARAM.
C
SUBROUTINE REGRSC (N,KBAND,RBAR,BETA)
INTEGER N !Read-only
INTEGER KBAND !Read-only
DOUBLE PRECISION RBAR(*) !Read-only
DOUBLE PRECISION BETA(*)
C
C JMH - 11/79 ANS FORTRAN 66
C
SUBROUTINE COVC (NCLMN,N,KBAND,RBAR,A)
INTEGER NCLMN !Read-only
INTEGER N !Read-only
INTEGER KBAND !Read-only
DOUBLE PRECISION RBAR(*) !Read-only
DOUBLE PRECISION A(NCLMN,NCLMN)
C
C JMH - 1/81 ANS FORTRAN 66
C
C IT - IT<0 IF ROW HAS ZEROES. IT>0 IF ROW HAS NO ZEROES.
C ILZ - INDEX OF LAST ROW WITH ZEROES.
C NRZ0 - TOTAL NUMBER OF ELEMENTS IN ROWS WITH ZEROES.
C IB - MAGIC NUMBER USED TO COMPUTE INDEX (NT) OF LAST
C ELEMENT IN ROW.
C ND - INDEX OF FIRST ELEMENT IN ROW.
C NT - INDEX OF LAST ELEMENT IN ROW.
C I - ROW INDEX.
C
C RBAR I ND NT IT
C
C D R R R R 0 0 0 0 0 0 T 1 1 6 -5 N = 10
C D R R R R 0 0 0 0 0 T 2 7 12 -4 KBAND = 4
C D R R R R 0 0 0 0 T 3 13 18 -3
C D R R R R 0 0 0 T 4 19 24 -2
C D R R R R 0 0 T 5 25 30 -1
C D R R R R 0 T 6 31 36 0 NRZ0 = 36
C D R R R R T 7 37 42 1
C D R R R T 8 43 47 2
C D R R T 9 48 51 3
C D R T 10 52 54 4
C
SUBROUTINE COV (NCLMN,N,A)
INTEGER NCLMN !Read-only
INTEGER N !Read-only
DOUBLE PRECISION A(NCLMN,NCLMN)
C
C COMPUTE UNSCALED COVARIANCE MATRIX ((A**T)*A)**(-1)
C
SUBROUTINE POLBAS (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(*)
C
C JMH - 9/88
C
C POLBAS COMPUTES THE FIRST N MONOMIALS AT X, EG. 1.0,X,X**2 ETC.
C THE RESULTS ARE RETURNED IN F. IX IS NOT USED. SEE FTL1L2 FOR AN
C EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION ROUTINES OF
C THIS FORM.
C
SUBROUTINE LEGBAS (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(*)
C
C JMH - 9/88
C
C LEGBAS COMPUTES THE FIRST N LEGENDRE POLYNOMIALS AT X, EG.
C 1.0, X, (3*X**2-1 )/2 ETC. THE RESULTS ARE RETURNED IN F.
C IX IS NOT USED BY LEGBAS. SEE FTL1L2 FOR AN
C EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION ROUTINES OF
C THIS FORM.
C THESE POLYNOMIALS ARE ORTHOGONAL ON THE INTERVAL [-1,+1]
C AND AS A RESULT ARE A BETTER CHOICE THAN THE MONOMIALS FOR
C LEAST SQUARES FITS TO DATA MORE OR LESS EVENLY DISTRIBUTED OVER
C THIS INTERVAL.
C
SUBROUTINE HRMBAS (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(*)
C
C JMH - 9/88
C
C HRMBAS COMPUTES THE FIRST N TERMS OF A HARMONIC EXPANSION AT X, EG
C 1.0, COS(X), SIN(X), COS(2*X), SIN(2*X), ... THE RESULTS ARE
C RETURNED IN F. IX IS NOT USED.
C
SUBROUTINE SPLBAS (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(NDIM,*)
C
C JMH - 9/88
C JMH - 7/02 - UPDATED
C
C SPLBAS COMPUTES THE N B-SPLINES OF ORDER K AND KNOT SEQUENCE K AT
C X. THE SPLINE AND ITS DERIVATIVES ARE RETURNED IN F. SEE FTL1L2
C FOR AN EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION
C ROUTINES OF THIS FORM. COMMON/SPFTCM/ (IN l1l2.h) MUST BE SET UP
C BEFORE SPLBAS IS CALLED. SEE SUBROUTINE KNOTS1. SPLBAS CALLS
C SUBROUTINES INTERV AND BSPLVD WHICH ARE FROM CARL DEBOOR'S
C B-SPLINE PACKAGE.
C
SUBROUTINE SPLBASP (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(NDIM,*)
C
C JMH - 2/02
C
C SPLBASP COMPUTES THE N PERIODIC B-SPLINES OF ORDER K AND KNOT
C SEQUENCE K AT X. THE RESULTS OR ONE OF THE FIRST THREE
C DERIVATIVES ARE RETURNED IN F. IX SPECIFIES THE ORDER OF THE
C DERIVATIVE, I.E IX=1 FOR THE FIRST DERIVATIVE. SEE FTL1L2 FOR AN
C EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION ROUTINES OF
C THIS FORM. COMMON/SPFTCM/ MUST BE SET UP BEFORE SPLBAS IS
C CALLED. SEE SUBROUTINE KNOTSP IN THE MATH LIBRARY FOR AN EXAMPLE
C OF A ROUTINE WHICH DOES THIS. SPLBASP CALLS SUBROUTINES INTERV
C AND BSPLVD WHICH ARE FROM CARL DEBOOR'S B-SPLINE PACKAGE AND ARE
C ALSO IN THE MATH LIBRARY.
C
SUBROUTINE SPLBASI (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(NDIM,*)
C
C JMH - 1/04
C
C SPLBASI COMPUTES THE N B-SPLINES OF ORDER K AND KNOT SEQUENCE K AT
C X. THE SPLINE, ITS DERIVATIVES AND A DEFINITE INTEGRAL DEFINED
C BY splbas.h ARE RETURNED IN F. SEE FTL1L2
C FOR AN EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION
C ROUTINES OF THIS FORM. COMMON/SPFTCM/ (IN l1l2.h) MUST BE SET UP
C BEFORE SPLBAS IS CALLED. SEE SUBROUTINE KNOTS1. SPLBAS CALLS
C SUBROUTINES INTERV AND BSPLVD WHICH ARE FROM CARL DEBOOR'S
C B-SPLINE PACKAGE.
C
C N - NUMBER OF B-SPLINES = NUMBER OF FIT COEFFICIENTS
C IX - INDEX OF THIS ELEMENT OF DATA ARRAY. NOT USED IN SPLBASI
C X - INDEPENDENT VARIABLE
C F(*,1) - VALUE OF B-SPLINE AT X
C F(*,2) - FIRST DERIVATIVE OF B-SPLINE AT X
C F(*,3) - SECOND DERIVATIVE OF B-SPLINE AT X
C F(*,4) - THIRD DERIVATIVE OF B-SPLINE AT X
C F(*,5) - INEGRAL OF B-SPLINE AROUND AT X AS DEFINED BY INTEGRATION
C COEFFICIENTS IN intdef.h
C
SUBROUTINE SPLBAS22 (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(NDIM,*)
C
C JMH - 7/02
C
C SPLBAS22 COMPUTES THE 2-DIMENSIONAL TENSOR PRODUCT OF BASIS
C FUNCTIONS SPLBAS AND SPLBAS.
C
SUBROUTINE SPLBASP2 (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(NDIM,*)
C
C JMH - 7/02
C
C SPLBASP2 COMPUTES THE 2-DIMENSIONAL TENSOR PRODUCT OF BASIS
C FUNCTIONS SPLBASP AND SPLBAS.
C
SUBROUTINE KNOTS1 (K,N,X1,X2,T)
INTEGER K !Read-only
INTEGER N !Read-only
DOUBLE PRECISION X1 !Read-only
DOUBLE PRECISION X2 !Read-only
DOUBLE PRECISION T(*)
C
C JMH - 9/88
C
C F:KNOTS1 - COMPUTES AN ARRAY OF KNOTS FOR A SPLINE FUNCTION.
C THE FIRST AND LAST INTER-KNOT SPACINGS ARE TWICE AS
C LONG AS THE EQUAL INTERIOR KNOT SPACINGS.
C
C INPUT PARAMETERS:
C K - ORDER (=DEGREE+1) OF B-SPLINES (K=4 FOR CUBIC SPLINES)
C X1 - LOCATION OF FIRST SPLINE KNOT
C X2 - LOCATION OF LAST SPLINE KNOT
C N - NUMBER OF B-SPLINES OF ORDER K FOR THE GIVEN KNOT
C SEQUENCE (N = NBKPT+K-2).
C
C OUTPUT PARAMETERS
C T - ARRAY OF KNOTS
C
C OTHER VARIABLES
C NBKPT - NUMBER OF BREAKPOINTS
C
C EXAMPLE (K=4, NBKPT=5, N=7):
C
C I1=1 I2=8
C | : | | | : |
C | : | <-DT-> | | : |
C | : | | | : |
C 1 5 6 7 8
C 2 9
C 3 10
C 4 11
C
SUBROUTINE KNOTS2 (K,N,X1,X2,T)
INTEGER K !Read-only
INTEGER N
DOUBLE PRECISION X1 !Read-only
DOUBLE PRECISION X2 !Read-only
DOUBLE PRECISION T(*)
C
C JMH - 9/88
C JMH - 11/01 MODIFIED FROM KNOTS
C
C F:KNOTS2 - COMPUTES AN ARRAY OF KNOTS FOR A SPLINE FUNCTION.
C THE FIRST AND LAST INTER-KNOT SPACINGS ARE HALF AS
C LONG AS THE EQUAL INTERIOR KNOT SPACINGS.
C
C INPUT PARAMETERS:
C K - ORDER (=DEGREE+1) OF B-SPLINES (K=4 FOR CUBIC SPLINES)
C X1 - LOCATION OF FIRST SPLINE KNOT
C X2 - LOCATION OF LAST SPLINE KNOT
C
C OUTPUT PARAMETERS
C T - ARRAY OF KNOTS
C N - NUMBER OF B-SPLINES OF ORDER K FOR THE GIVEN KNOT
C SEQUENCE (N = NBKPT+K-2).
C
C OTHER VARIABLES
C NBKPT - NUMBER OF BREAKPOINTS
C
C EXAMPLE (K=4, NBKPT=5, N=7):
C
C I1=1 I2=8
C | | | | |
C | | <-DT-> | | |
C | | | | |
C 1 5 6 7 8
C 2 9
C 3 10
C 4 11
C
SUBROUTINE KNOTSP (K,N,X1,X2,T)
INTEGER K !Read-only
INTEGER N !Read-only
DOUBLE PRECISION X1 !Read-only
DOUBLE PRECISION X2 !Read-only
DOUBLE PRECISION T(*)
C
C JMH - 6/02
C
C F:KNOTSP - COMPUTES AN ARRAY OF KNOTS FOR A PERIODIC SPLINE
C FUNCTION.
C
C INPUT PARAMETERS:
C K - ORDER (=DEGREE+1) OF B-SPLINES (K=4 FOR CUBIC SPLINES)
C X1 - LOCATION OF FIRST SPLINE KNOT
C X2 - LOCATION OF LAST SPLINE KNOT
C
C OUTPUT PARAMETERS
C T - ARRAY OF KNOTS
C N - NUMBER OF B-SPLINES OF ORDER K FOR THE GIVEN KNOT
C SEQUENCE (N = NBKPT+K-2).
C
C OTHER VARIABLES
C NBKPT - NUMBER OF BREAKPOINTS
C
C EXAMPLE (K=4, NBKPT=5, N=7):
C
C I1=1 I2=8
C | | | | |
C | | <-DT-> | | |
C | | | | |
C 1 5 6 7 8
C 2 9
C 3 10
C 4 11
C
SUBROUTINE INTERV (XT,LXT,X,ILEFT,MFLAG)
DOUBLE PRECISION XT(LXT) !Read-only
INTEGER LXT !Read-only
DOUBLE PRECISION X !Read-only
INTEGER ILEFT
INTEGER MFLAG
C OMPUTES LARGEST ILEFT IN (1,LXT) SUCH THAT XT(ILEFT) .LE. X
SUBROUTINE BSPLVD (T,K,X,ILEFT,VNIKX,NDERIV)
DOUBLE PRECISION T(*) !Read-only
INTEGER K !Read-only
DOUBLE PRECISION X !Read-only
INTEGER ILEFT !Read-only
DOUBLE PRECISION VNIKX(K,NDERIV)
INTEGER NDERIV !Read-only
C CALCULATES VALUE AND DERIV.S OF ALL B-SPLINES WHICH DO NOT VANISH
C AT X. FILL VNIKX(J,IDERIV), J=IDERIV, ... ,K WITH NONZERO VALUES
C OF B-SPLINES OF ORDER K+1-IDERIV , IDERIV=NDERIV, ... ,1, BY
C REPEATED CALLS TO BSPLVN
C DIMENSION T(ILEFT+K)
SUBROUTINE BSPLVN (T,JHIGH,INDEX,X,ILEFT,VNIKX)
DOUBLE PRECISION T(*) !Read-only
INTEGER JHIGH !Read-only
INTEGER INDEX !Read-only
DOUBLE PRECISION X !Read-only
INTEGER ILEFT !Read-only
DOUBLE PRECISION VNIKX(JHIGH)
C ALCULATES THE VALUE OF ALL POSSIBLY NONZERO B-SPLINES AT 'X' OF
C ORDER MAX(JHIGH,(J+1)(INDEX-1)) ON 'T'.
C DIMENSION T(ILEFT+JHIGH)
SUBROUTINE LPBASIS (N,IX,X,F)
INTEGER N !Read-only
INTEGER IX !Read-only
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION F(NDIM,5)
C
C JMH - 1/04
C
C LPBASIS COMPUTES THE N B-SPLINES OF ORDER K AND KNOT SEQUENCE K AT
C X. THE SPLINE, ITS DERIVATIVES AND A DEFINITE INTEGRAL DEFINED
C BY splbas.h ARE RETURNED IN F. SEE FTL1L2
C FOR AN EXAMPLE OF A SUBROUTINE WHICH CALLS BASIS FUNCTION
C ROUTINES OF THIS FORM. COMMON/SPFTCM/ (IN l1l2.h) MUST BE SET UP
C BEFORE SPLBAS IS CALLED. SEE SUBROUTINE KNOTS1. LPBASIS CALLS
C SUBROUTINES INTERV AND BSPLVD WHICH ARE FROM CARL DEBOOR'S
C B-SPLINE PACKAGE.
C
C N - NUMBER OF B-SPLINES = NUMBER OF FIT COEFFICIENTS
C IX - INDEX OF THIS ELEMENT OF DATA ARRAY. NOT USED IN LPBASIS
C X - INDEPENDENT VARIABLE
C F(*,1) - VALUE OF B-SPLINE AT X
C F(*,2) - FIRST DERIVATIVE OF B-SPLINE AT X
C F(*,3) - SECOND DERIVATIVE OF B-SPLINE AT X
C F(*,4) - THIRD DERIVATIVE OF B-SPLINE AT X
C F(*,5) - INEGRAL OF B-SPLINE AROUND AT X AS DEFINED BY
C INTEGRATION COEFFICIENTS IN intdef.h
C
SUBROUTINE MAKENOTES ()
C
SUBROUTINE MINV (A,N,D,L,M)
DOUBLE PRECISION A(*)
INTEGER N !Read-only
DOUBLE PRECISION D
INTEGER L(*)
INTEGER M(*)
C
C ...............................................................
C
C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C STATEMENT WHICH FOLLOWS.
C
C DOUBLE PRECISION A,D,BIGA,HOLD
C
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C ROUTINE.
C
C THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS. ABS IN STATEMENT
C 10 MUST BE CHANGED TO DABS.
C
C ...............................................................
C
C SEARCH FOR LARGEST ELEMENT
C
SUBROUTINE OFAM2ATM ()
C
C JMH - 10/93
C
SUBROUTINE OFAM2CEDAR ()
C
C OFAM2CEDAR TRANSFERS THE CONTENTS OF THE INSCAL OUTPUT STRUCTURE,
C OFAM, TO A CEDAR RECORD AND WRITES THE RECORD TO A CEDAR FILE.
C
SUBROUTINE PARAM (ASET,IACF,A)
DOUBLE PRECISION ASET(*) !Read-only
INTEGER IACF !Read-only
DOUBLE PRECISION A(*)
C
C JMH - 9/81
C JMH - 10/83 COPIED FROM 2INSC*F:PARAM. THIS WILL BE USED AS
C STARTING POINT FOR USING NEW DATABASE ROUTINES IN INSCAL.
C
C PARAM CALCULATES INITIAL GUESSES FOR THE ACF PARAMETERS A.
C IF ASET CONTAINS A VALUE OTHER THAN RMISS, THAT VALUE IS
C USED. OTHERWISE, PARMOD IS CALLED.
C
C ASET - SPECIFIED VALUE FOR A PARAMETER
C IACF - INDEX OF ACF
C A - INITIAL GUESS FOR PARAMETER
C
DOUBLE PRECISION FUNCTION PARMOD (IACF,I)
INTEGER IACF !Read-only
INTEGER I !Read-only
C
C JMH - 9/81
C JMH - 10/83 COPIED FROM 2INSC*F:PARMOD
C GBL - 5/84 CHANGED CALC. OF "A" FOR E-DENSITY
C
SUBROUTINE A2SP (IT,A1,A2,SP1,SP2)
INTEGER IT !Read-only
DOUBLE PRECISION A1 !Read-only
DOUBLE PRECISION A2 !Read-only
DOUBLE PRECISION SP1
DOUBLE PRECISION SP2
C
SUBROUTINE SP2A (IT,SP1,SP2,A1,A2)
INTEGER IT !Read-only
DOUBLE PRECISION SP1 !Read-only
DOUBLE PRECISION SP2 !Read-only
DOUBLE PRECISION A1
DOUBLE PRECISION A2
C
SUBROUTINE TITE2SP (IT,TI,TE,SP1,SP2)
INTEGER IT !Read-only
DOUBLE PRECISION TI !Read-only
DOUBLE PRECISION TE !Read-only
DOUBLE PRECISION SP1
DOUBLE PRECISION SP2
C
SUBROUTINE SP2TITE (IT,SP1,SP2,TI,TE)
INTEGER IT !Read-only
DOUBLE PRECISION SP1 !Read-only
DOUBLE PRECISION SP2 !Read-only
DOUBLE PRECISION TI
DOUBLE PRECISION TE
C
SUBROUTINE DSP2DT (IT,SP1,SP2,DSP1,DSP2,DSP12,DTI,DTE,DTR)
INTEGER IT !Read-only
DOUBLE PRECISION SP1 !Read-only
DOUBLE PRECISION SP2 !Read-only
DOUBLE PRECISION DSP1 !Read-only
DOUBLE PRECISION DSP2 !Read-only
DOUBLE PRECISION DSP12 !Read-only
DOUBLE PRECISION DTI
DOUBLE PRECISION DTE
DOUBLE PRECISION DTR
C
SUBROUTINE DSP2DA (IT,SP1,SP2,DSP1,DSP2,DSP12,DA1,DA2)
INTEGER IT !Read-only
DOUBLE PRECISION SP1 !Read-only
DOUBLE PRECISION SP2 !Read-only
DOUBLE PRECISION DSP1 !Read-only
DOUBLE PRECISION DSP2 !Read-only
DOUBLE PRECISION DSP12 !Read-only
DOUBLE PRECISION DA1
DOUBLE PRECISION DA2
C
SUBROUTINE A2TITE (A1,A2,TI,TE)
DOUBLE PRECISION A1 !Read-only
DOUBLE PRECISION A2 !Read-only
DOUBLE PRECISION TI
DOUBLE PRECISION TE
C
SUBROUTINE TITE2A (TI,TE,A1,A2)
DOUBLE PRECISION TI !Read-only
DOUBLE PRECISION TE !Read-only
DOUBLE PRECISION A1
DOUBLE PRECISION A2
C
SUBROUTINE POSTINT (QPOST,ISTAT)
INTEGER QPOST !Read-only
INTEGER ISTAT
C
C MEDIAN FILTERS SUB-INTEGRATION-PERIOD ACFS IN TIME
C
SUBROUTINE PRECSM (LFN)
INTEGER LFN !Unknown (program incomplete)
C
C JMH - 12/92
C
C PRECSM WRITE A ONE-LINE SUMMARY OF THE CURRENT DATA RECORD
C TO THE SPECIFIED LOGICAL UNIT.
C
SUBROUTINE PROFINT (OASISDIR,ISTAT)
CHARACTER OASISDIR !Unknown (program incomplete)
INTEGER ISTAT
C
C PROFINT IMPLEMENTS ALTITUDE INTEGRATION OF THE IFAM ACF ARRAY.
C BY DEFAULT, THE ACF ARRAY IS NOT MODIFIED. TO SPECIFY ALTITUDE
C INTEGRATION FOR A SPECIFIED MODE, INCLUDE A FILE OF THE FORM
C MODENAME.DAT IN THE OASISDIR IN WHICH INSCAL IS RUN, E.G.
C A480.30M8910.DAT. THIS FILE HAS THE FOLLOWING FORMAT:
C
C MDNAME = A480.30M8910
C MDPTYPE = 2
C MDNROWS = 5
C MDNCOLS = 3
C 150.0 200.0 3.0
C 200.0 300.0 6.0
C 300.0 400.0 8.0
C 400.0 500.0 10.0
C 500.0 800.0 12.0
C MDPTYPE = 1 (N-POINT CENTERED MEAN), MDTYPE = 2 (N-POINT
C CENTERED MEDIAN), MDTYPE = 3 (L2 SPLINE FIT) AMD MDPTYPE= 4 (L1
C SPLINE FIT) ARE SUPPORTED. THE INTEGRATION PARAMETERS ARE
C MDPTYPE DEPENDENT AND ARE ENTERED AS AN ARRAY WITH MDROWS ROWS
C AND MDCOLS COLUMNS.
C
C FOR MDPTYPE = 1 AND 2, EACH ROW CONTAINS A MINIMUM AND MAXIMUM
C ALTITUDE AND THE A SPECIFICATION OF THE NUMBER OF POINTS IN THE
C MEAN OR MEDIAN. THUS THE FIRST ROW STATES THAT BETWEEN 200 AND
C 300 KM A 13-POINT AVERAGE WILL BE COMPUTED. ACFS BELOW 150 KM
C ARE LEFT UNCHANGED IN THIS EXAMPLE.
C
C FOR MDTYPE = 3 AND 4, THE FORMAT IS AS FOLLOWS
C
C MDNAME = A480.30M8910
C MDPTYPE = 4
C MDNROWS = 1
C MDNCOLS = 4
C 100.0 5.0 600.0 100.0
C
C IN THE EXAMPLE,
C
C Z1 = 100.0 - LOWEST SPLINE KNOT
C DZ1 = 5.0 - DISTANCE FROM Z1 TO NEXT SPLINE KNOT
C Z2 = 600.0 - HIGHEST SPLINE KNOT
C DZ2 = 100.0 - APPROXIMATE KNOT SPACING AT Z2
C
C IN PRACTICE THE ACTUAL SPACING BETWEEN THE SECOND-HIGHEST KNOT
C AND Z2 WILL BE SOMEWHAT LARGER THAN DZ2. KNOT SPACING INCREASES
C LINEARLY WITH ALTITUDE. THE EXAMPLE YIELDS THE FOLLOWING KNOTS:
C
C 100.0 105.0 111.2 118.6 127.5 138.1 150.9 166.2 184.5
C 206.4 232.6 264.0 301.7 346.8 400.8 465.5 600.0
C
SUBROUTINE PRPARM (LFNLIS)
INTEGER LFNLIS !Read-only
C
C JMH - 8/99
C
C PRPARM PRINTS THE MHR PARAMETER FILE.
C
SUBROUTINE PRSFIL (OUTFILE,LFN)
CHARACTER OUTFILE !Unknown (program incomplete)
INTEGER LFN !Unknown (program incomplete)
C
C JMH - 10/00
C
C PRECSM WRITE A ONE-LINE SUMMARY OF THE CURRENT DATA RECORD
C TO THE SPECIFIED LOGICAL UNIT. THIS FUNCTION IS USED TO
C WRITE A ONE LINE STATUS INSCAL STATUS FILE WHICH CONTAINS
C INFORMATION ABOUT THE LAST RECORD WRITTEN BY INSCAL AND THE
C UNIX TIME WHEN IT WAS LAST WRITTEN.
C
SUBROUTINE SETFIT (NLAGI,NLAGP2I,NCORI,NDCMTI,NLAGSCI,IFTFNI,DTAUI,ACFCRI)
INTEGER NLAGI !Read-only
INTEGER NLAGP2I !Read-only
INTEGER NCORI !Read-only
INTEGER NDCMTI !Read-only
INTEGER NLAGSCI !Read-only
INTEGER IFTFNI !Read-only
DOUBLE PRECISION DTAUI !Read-only
DOUBLE PRECISION ACFCRI(NCORMX,NLAGMX) !Read-only
C
C SETFIT SETS THE VARIABLES IN FITC. THESE VARIABLES ARE NEEDED BY
C FUNCTN, WHICH COMPUTES THE THEORETICAL ACF.
C
SUBROUTINE SETCON (IER)
INTEGER IER
C
C JMH - 9/81
C
C SETCON EVALUATES THE CONSTANTS REQUIRED IN THE CALCULATION OF AN
C ARBITRARY INCOHERENT SCATTER AUTOCORRELATION FUNCTION. ALL
C CONSTANTS ARE IN MKSA UNITS. FREQ, B AND ALPHA MUST BE STORED
C IN COMMON/DATA/ BEFORE SETCON IS CALLED. THE REMAINING
C CONSTANTS ARE CALCULATED BY SETCON AND STORED IN DATA.
C
C FREQ - RADAR FREQUENCY (MHZ)
C B - MAGNETIC FIELD (WEBER M-2)
C ALPHA - ANGLE BETWEEN WAVE VECTOR AND B (RAD)
C RE - ELECTRON RADIUS (M)
C AMU - ATOMIC MASS UNIT (KG)
C E - ELECTRON CHARGE (C)
C BOLTK - BOLTZMAN CONSTANT (J DEG-1)
C RME - ELECTRON MASS (KG)
C PI - 3.1415926536
C RLAM - RADAR WAVELENGTH (M)
C
C THE REMAINING FUNCTIONS EVALUATED BY SETCON ARE FUNCTIONS OF
C THESE CONSTANTS.
C
C IER = 201 FREQ .LE. 0. FREQ SET EQUAL TO 440.2 MHZ.
C
SUBROUTINE SETSTAT (SN,QGOOD,NIT,IF)
DOUBLE PRECISION SN !Read-only
LOGICAL QGOOD !Read-only
INTEGER NIT !Read-only
INTEGER IF !Read-only
C
SUBROUTINE SPECT2 (TI,TE,DF,P2,RM1,RM2,INDEL,N,DELOM,SP)
DOUBLE PRECISION TI !Read-only
DOUBLE PRECISION TE !Read-only
DOUBLE PRECISION DF !Read-only
DOUBLE PRECISION P2 !Read-only
DOUBLE PRECISION RM1 !Read-only
DOUBLE PRECISION RM2 !Read-only
INTEGER INDEL !Read-only
INTEGER N !Read-only
DOUBLE PRECISION DELOM !Read-only
DOUBLE PRECISION SP(*)
C
C JMH - 9/81
C
C SPECT2 CALCULATES THE COLLISIONLESS INCOHERENT SCATTER SPECTRUM
C CORRESPONDING TO A BINARY MIXTURE OF IONS. THIS VERSION IS
C ACCURATE TO 11 SIGNIFICANT FIGURES.
C TI - ION TEMPERATURE
C TE - ELECTRON TEMPERATURE
C DF - DEBYE FACTOR
C P2 - PERCENTAGE OF SECOND ION
C RM1 - MASS OF FIRST ION
C RM2 - MASS OF SECOND ION
C INDEL - ELECTRON ADMITTANCE CALCULATED IF INDEL .GT. 2
C N - NUMBER OF SPECTRUM POINTS TO BE CALCULATED
C OMEGA - ANGULAR FREQUENCIES OF SPECTRUM POINTS
C SP - SPECTRUM
C
SUBROUTINE SPECTC (T,CL,N,DELOM,SP)
DOUBLE PRECISION T !Read-only
DOUBLE PRECISION CL !Read-only
INTEGER N !Read-only
DOUBLE PRECISION DELOM !Read-only
DOUBLE PRECISION SP(*)
C
C JMH - 9/84 HARRIS FORTRAN 77
C 3/85 MODIFIED TO USE DELOM INSTEAD OF OMEGA AND TO OUTPUT
C SPECTRUM NORMALIZED TO 1 AT ZERO LAG.
C
C SPECTC CALCULATES THE INCOHERENT SCATTER SPECTRUM WITH COLLISIONS
C AS DESCRIBED IN DOUGHERTY AND FARLEY, JGR, 68, 5473-5486, 1963.
C
SUBROUTINE SPLITMF (MFFULL,MFDIR,MF)
CHARACTER MFFULL !Unknown (program incomplete)
CHARACTER MFDIR
CHARACTER MF
C
SUBROUTINE STACOR (KINST,SLATGD,SLON,SALTGD,SLATGC,SR,ISTAT)
INTEGER KINST !Read-only
DOUBLE PRECISION SLATGD
DOUBLE PRECISION SLON
DOUBLE PRECISION SALTGD
DOUBLE PRECISION SLATGC
DOUBLE PRECISION SR
INTEGER ISTAT
C
C Returns Instrument location coordinates from instrument code.
C
C Input:
C KINST - Instrument code.
C
C Output:
C SLATGD - Instrument geodetic latitude
C SLON - Instrument longitude
C SALTGD - Instrument geodetic latitude
C SLATGC - Instrument geocentric latitude
C SR - Radial distance of instrument from center of earth
C ISTAT - 0 if successful.
C
DOUBLE PRECISION FUNCTION VELX (ACFRR,ACFII,NLAGS,DELT,FREQ)
DOUBLE PRECISION ACFRR(*) !Read-only
DOUBLE PRECISION ACFII(*) !Read-only
INTEGER NLAGS !Read-only
DOUBLE PRECISION DELT !Read-only
DOUBLE PRECISION FREQ !Read-only
C
SUBROUTINE WOFZ (X,Y,WRE,WIM)
DOUBLE PRECISION X !Read-only
DOUBLE PRECISION Y !Read-only
DOUBLE PRECISION WRE
DOUBLE PRECISION WIM
C
C ACM ALGORITHM 363
C WALTER GAUTSCHI
C FORTRAN TRANSLATION - J.M.HOLT - 9/15/76
C 9/17/84 - SOLUTION EXTENDED TO WHOLE Z-PLANE.
C
C THIS PROCEDURE EVALUATES THE REAL AND IMAGINARY PART OF THE
C FUNCTION W(Z) = EXP(-Z**2)*ERFC(-I*Z) FOR ARGUMENTS Z = X+I*Y IN
C THE FIRST QUADRANT OF THE COMPLEX PLANE. THE ACCURACY IS 10
C DECIMAL PLACES AFTER THE DECIMAL POINT, OR BETTER. FOR THE
C UNDERLYING ANALYSIS, SEE W. GAUTSCHI, 'EFFICIENT COMPUTATION OF
C
Notes
+ Data structures in inscal
Inscal global parameters are defined in inscal.h
Inscal uses common blocks to organize data. These common blocks are
defined in include files. The key common blocks are:
COMMON /IFAM/ - ifamf.h
acfdprn.f
acfit.f
ckifam.f
efitfp.f
fitstats.f
fixifam.f
gtiftm.f
if2ofi.f
ifam_in_f.f
ifam_out_f.f
inscal.f
makehtml.f
parmod.f
precsm.f
prsfil.f
COMMON /OFAM/ - ofamf.h
acfit.f
atm2ofam.f
fitstats.f
if2ofi.f
inscal.f
makehtml.f
ofam2atm.f
ofam2cedar.f
ofam_in_f.f
ofam_out_f.f
plotall.f
plotne.f
plotpar.f
plottite.f'
plotvo.f
COMMON /INSCOM/ - inscom.h
exinsc.f
inscal.f
precsm.f
prsfil.f
COMMON /CNTRL/ - cntrl.h
cmdin.f
cmdout.f
exinsc.f
gtcmtm.f
inicom.f
inscal.f
makenotes.f
COMMON /SSWCOM/ - qssweq.h
acfit.f
exinsc.f
fixifam.f
inicom.f
inscal.f
parmod.f
COMMON /PARCOM/ - qssweq.h
acfit.f
exinsc.f
fixifam.f
inicom.f
inscal.f
parmod.f
COMMON /FITC/ - fitc.h
acfit.f
fitg.f
functn.f
COMMON /PARTAB/ - partab.h
acfit.f
param.f
parini.f
COMMON /FTSTAT/ - stat.h
acfit.f
inicom.f
prcmst.f
prstat.f
setstat.f
Character data are in a separate Common block in accordance with the
ANSII Fortran standard. This ifamf.h has COMMON /IFAMC/ in addition to
COMMON /IFAM/. Common blocks also contain save statements, e.g.
SAVE /GEOCOM/ to ensure that Common block variables are retained
when the Common block is not referenced by the main program.
+ ACFs in inscal
Inscal gets ACFs from the ifam structure, which is defined in ifamf.h
and ifamc.h. Ifam contains everything inscal knows about an I. S.
dataset. Note in particular that inscal only deals with ACFs, not the
more general lag-products.
ACFs contain a maximum of NACFMX lags. The actual number is given by
ifam array variable NLAGS, which contains separate values for each
ACF. Acfit uses array index variable IACF for the ACFs in ifam, and
most of the code in acfit is inside a "DO IACF=1,NACF" loop. Inscal
deals with instrumental corrections by modifying theoretical ACF's to
match measured ACFs. It is usually desirable to use a finer lag
spacing for these corrections than present in the measured ACF. In this
case, the correction matrix used to correct theoretical ACFs contains a
factor of NDCMT more lags than the measured ACF. The number of lags in
the correction matrix is given by ifam array variable NLAGSC, which
contains separate values for each ACF. After the correction have been
calculated, the theoretical ACF is decimated by a factor of NDCMT.
The fitting program requires the theoretical ACFs to have a number of
lags which is equal to a power of 2. Acfit sets the variable NLAGP2 to
the first power of two greater than or equal to NLAGSC(IACF), and all
teoretical ACFs are calculated with NLAGP2 lags. NLAGP2 is passed to
functn, which calculates the theoretical ACF, in fitc.h. The measured
ACF is zero-padded from NLAGS to NLAGP2/NDCMT. The measured ACF from
ifam is transfered to a local variable ACF in acfit, which passes it as
an argument to fit1, which passes it as an argument to curfit, which
passes it as an argument to functn. ACF contains both the real and
imaginary parts of the ACF. The real part starts at index 1 and the
imaginary part starts at index NLAG2+2, where NLAG2 = NLAGP2/NDCMT1.
Acfit also defines the scalar variables NLAG=NLAGS(IACF),
NLAGSC1=NLAGSC(IACF) and M = NLAG2+1.
+ External Functions and Subroutines called by inscal
Inscal requires the Madrec, Geophysics and NCAR Graphics libraries
The Fortran Library calls may create portability problems.
AGSETC - NCAR Graphics
AGSETF - NCAR Graphics
AGSETI - NCAR Graphics
AGSETP - NCAR Graphics
CEDARCLOSE - Madrec
CEDARCREATE - Madrec
CEDARGETSCALEFACTOR - Madrec
CEDAROPEN - Madrec
CEDARSET1DPARM - Madrec
CEDARSET2DPARM - Madrec
CEDARWRITE - Madrec
CLOSE_MIDAS_FILE - midas_io.c
CONVRT - Geophysics
DATER - Geophysics
DTIME - Fortran
EXIT - Fortran
EZMXY - NCAR Graphics
FDATE - Fortran
FLUSH - Fortran
FRAME - NCAR Graphics
GACWK - NCAR Graphics
GCLKS - NCAR Graphics
GCLWK - NCAR Graphics
GDAWK - NCAR Graphics
GETARG - Fortran
GETCWD - Fortran
GETPID - Fortran
GOPKS - NCAR Graphics
GOPWK - NCAR Graphics
GPM - NCAR Graphics
GSCHH - NCAR Graphics
GSCHSP - NCAR Graphics
GSCHXP - NCAR Graphics
GSMK - NCAR Graphics
GSMKSC - NCAR Graphics
GSPLCI - NCAR Graphics
GSTXAL - NCAR Graphics
GSTXCI - NCAR Graphics
GTS5 - Geophysics (MSIS)
GTX - NCAR Graphics
HOSTNM - Fortran
IARGC - Fortran
IEEE_FLAGS - Fortran
JDATER - Geophysics
JDAY - Geophysics
LNBLNK - Fortran
METERS - Geophysics (MSIS)
NGSETC - NCAR Graphics
OPEN_MHR_PARMS - midas_io.c
OPEN_MIDAS_FILE - midas_io.c
PLOTIF - NCAR Graphics
POINT - Geophysics
READ_MIDAS_RECORD - midas_io.c
SATELLITE_DETECT - satellite_detect.c
SET - NCAR Graphics
SIGNAL - Fortran
TIME - Fortran
TRETRV - Geophysics (MSIS)
TSELEC - Geophysics (MSIS)
TO DO
_____
Verify how ACF errors, variances are used (DELACF, ACFVAR).