* Tue Jun 21 14:11:02 EDT 1994 /plas1/h1/mzh/src/test mjsanl.f
C Mon Oct 18 14:45:37 EDT 1993 /plas5/h1/vgr/src/mjs mjsanl
* Fri Jan 11 15:22:10 EST 1991 /herschel/h1/rlm/src/test mjsanl.f
c
c XFULL has been modified as well as MJSANL to use
c full analysis on GS5 Neptune data - checkout is
c still in progress (as well as debugging!)
c
c 10:30 AM 4/23/90 rlm
c
c
c Begin final modifications and check out to support
c Neptune encounter operations
c RLM 3:25 PM Saturday July 15, 1989
c
c Modifications to support thermal cycle
c testing on Voyager 1
c
c January 5, 1989
c
c Begin making modifications to use PMODEL to handle eletron to ion
c and ion to electron feedthrough in fits
c 2:12 PM, Sunday, October 2, 1988 - rlm
c
c
COMMON/PLTSIZ/ W,H
COMMON/ANSWER/ANS(150),CURRNT(512)
common/full/curnew(544),ifull,jnei
common/sedr/isedr,ist(6),isp(6)
COMMON/SPEC/SBLK(60)
INTEGER*2 JTB(6),JDATA(512),JNE,JTLMOD,JCLK
DIMENSION TEMP(3)
integer i
integer isedr
character*4 itype
call initio( 1, 36,i)
call initio( 1,-13,i)
call initio( 1, 29,i)
3 format(A130)
4 read (29,3,end=5)
go to 4
5 continue
c
c Initialize numerical data output file --- asz 10/25/93
c
call initio( 1, 28,i)
6 read (28,3,end=7)
go to 6
7 continue
c
CALL NXTMOD (ANS,JTB,JDATA,JNE,LSTAT,JTLMOD,JCLK,TEMP,IEOD,IL,
*ITYPE)
c
call mjsanl(jtb)
end
SUBROUTINE MJSANL(jtb)
c
c ********************************************************************
c * *
c * MAIN *
c * *
c * GENERAL MAGNETOSPHERIC PLASMA DATA ANALYSIS PROGRAM: *
c * *
c * AN INTERACTIVE SYSTEM FOR PERFORMING NON-LINEAR LEAST *
c * SQUARES FITS TO DATA OBTAINED BY THE "PLASMA SUBSYSTEM" *
c * (PLS) EXPERIMENT ON BOARD THE VOYAGER 1 AND 2 SPACECRAFT *
c * DURING THEIR ENCOUNTERS WITH THE OUTER PLANETS JUPITER, *
c * SATURN, URANUS, AND NEPTUNE. *
c * *
c * *
c * ALSO CONTAINS COMMON BLOCK "BMODEL" *
c * *
c ********************************************************************
c
DIMENSION ITB(6),KTB(6)
INTEGER*2 JTB(6),JDATA(512),JNE,JTLMOD,JCLK
c
c Differentiate between ntype and nntype to decide
c whether to read in initialization parameters from
c answer array (i.e., Mission Tape, or Summary or LONG
c spectral tape during cruise = use of PHASE3).
c
c June 13, 1990 RLM
c
integer*4 ntype
integer*4 nntype
c
c Input data type common block
c
common/tptype/jtype,ntype
c
c JTYPE C4: The data tape/file type string returned by NXTMOD
c in the last variable ITYPE. Legal possibilities are:
c 'EDR ','SUM ','SPL ', and 'VGRT' corresponding to
c EDR, summary, spectral, and mission tapes, respectively.
c NTYPE I4: Index number of the tape/file type: 1=EDR; 2=SUM;
c 3=SPL; 4=VGRT. Half way through mjsanl.f, it is reset
c to 4 for SUM and Long SPL files. NTYPE=4 then indicates
c data tapes/files with available ANS array.
c
c Add to use SETJTL 3/23/90
c
integer*2 jtlmo,jtloff
c
c 7/19/89 NET for Neptune
c
integer*2 njtb(6)
c
character*4 ifrmt1
COMMON/HEADER/PROJID,FILEID,ISCID,TAPEID(2),IGENDT,IGENTM,FIPID(2)
1,IFGEND,IFGENT,DPTJ1(2),IFRMT1(2),DPTJ2(2),IFRMT2(2),DPTJ3(2),IFRM
2T3(2),DPTJ4(2),IFRMT4(2),HSPARE(18)
c
c The only element of this common block used is the S/C id
c stored in ISCID. 0 = Voyager 2; 1 = Voyager 1. It is loaded
c in phase2 and is needed only in encounter mode.
c
integer iplt,i5,icur,irec,nrep
common /cntrl/ iplt,i5,icur,irec,nrep
c
c iplt Plotting flag (1=Yes, 0=No)
c i5 input file number, default 5, alt = 15
c icur number for type of analysis,
c if negative calculate currents
c irec record results on unit # 17
c nrep number of iterations befor ploting
c
COMMON/ANSWER/ANS(150),CURRNT(512)
c
c The ANS array is primarily filled in NXTMOD. Only the following
c elements are filled here:
c ANS(20): In encounter mode: The maximum response of any cup
c to rigid corotating plasma.
c In cruise mode: The proton number density copied from
c ANS(38).
c ANS(21): The cup number which gave the maximum response above.
c CURRNT(512): The real currents (fA) in each channel for each cup.
c The currents are calculated in KNTCUR and passed back.
c
common/crt/icrt,ns
c
c ICRT: Flag to correct the L mode currents for noise. 1=Yes, 0=No
c NS: The number of spectra to average. 0=no averaging
c
common/sedr/isedr,ist(6),isp(6)
c
c ISEDR: SEDR tape/file type index: 0=No SEDR available;
c 2=SSEDR (via phase2); 3=MSSEDR (via phase3)
c IST(6): The beginning of the SSEDR file (integer array)
c ISP(6): The end of the SSEDR file.
c
common/sel/isel
c
c ISEL: Index of the main mode of operation:
c 0 = Cruise mode; 1 = Encounter mode
c
COMMON/BMODEL/IBMODL,IDSRNN,IK(10)
c
c BMODEL is only initialized in mjsanl.
c
c Array curnew is for data in full up display (ions and electrons)
c
common/full/curnew(544),ifull,jnei
c
c CURNEW(544): real currents (fA) in each channel for each cup
c IFULL: 0=no joint ion/electron analysis; 1= Joint L,E1,E2 analysis;
c 2= Joint M,E1,E2 analysis
c JNEI: The # of channels in a cup for the current mode
c 128 for M, 16 for L,E1, and E2 modes
c
c Pass down step information for PMODEL
c
common/pstep/rstl,rstm,rste1,rste2,
+rstle,rstme,rste1i,rste2i
COMMON/MID/NUMFOR
c
c NUMFOR: The planet id: 0=Cruise mode; 1=Jupiter; 2=Saturn
c 3=Uranus; 4=Neptune
c
COMMON/PLTSIZ/W,H
c
c W: pagelength in inches
c H: pagewidth in inches
c
c Character strings to hold the name of the mission data
c and short sedr file names
c
character*64 ftdata,ftsedr
c
c Character strings to hold the start and stop times of the
c attached vgr and mssedr tapes
c
character*24 sedrtime, vgrtime
c
c Character string to hold type of run from the RUN_TYPE environment
c
character*3 run_type
c
dimension bcupn(4)
c
c Set up array for averaging currents
c
dimension curave(512)
dimension ncount(512)
c
c Common block for elapsed CPU time. Not used currently.
c
COMMON/TIMExx/ITT,IAA
COMMON/RTIM/RCUR(6)
c
c RCUR(6): The current six element time array of the
c most recently read spectra.
c
character*4 TIMESR(5)
COMMON/CDATES/SPARA(17),TIMESR
c
c There are many time related stuff in this common block,
c but mjsanl uses only TIMESR
c
c TIMESR(5): When printing side by side the 5 4-character
c strings they give the time of the current spectra
c in the format: '19xx xxx xxxx:xx.xxx'
c
COMMON/SCID/VOY,IME
c
c VOY: 453=Voyager 1; 372=Voyager 2
c IME: The number of channels in a cup for current mode.
c 128 for M, 16 for L,E1, and E2 modes
c
COMMON/SCALE/YB,IDEC,XB,XSCALE
c
c The coordinate limits of the current vs. channel # for each cup.
c YB: The Y axis lower limit (fA)
c IDEC: The number of logarithmic decades on the Y axis
c XB: The X axis lower limit in channel number
c XSCALE: The length of the X axis scale in channel numbers
c
COMMON/SPEC/SBLK(60)
c
c SBLK(60): Trajectory and coordinate transformation data
c filled in xcru(phase3.f) and xspn(phase2.f).
c
COMMON/PLCONS/SPARE(726),CN(4,3),COSB(3),COS2B(3),SIN2B(3)
* ,CALCUR(12),CLOCK(3),CAP(2),GAIN(2),V0,ALPHB
* ,NQ(4),MASS(4)
c
c Out of all of these general plasma analysis parameters,
c mjsanl uses only V0: 453.=Voy. 1; 372.=Voy.2
c
dimension cupn(3,4)
character*4 itype,jtype
character*4 idtype(4)
c
c Offset correction currents for L long (GS5) at Neptune
c
dimension curoff(64)
c
c SET UP CONSTANTS USED IN XSIMUL AND PMODEL
c
LOGICAL*4 LD
COMMON/NLLS/ICPLT,MP,ME,RSTEP,R1,R2,IPRT,LD
c
c ICPLT: Fitting mode id: 0=no fit; 1=PMODEL; 3=XJC; 4=ABSO;
c 5=ABSOM; 6=BiMax. ABSO; 7=Cold beam; 9=Very cold beam;
c >10=Terminate.
c MP: Spectral type: 0=L; 1=M; 2=E1; 3=E2
c ME: Spectral type: 0=E2; 1=E1; 2=M; 3=L
c RSTEP,R1,R2: PMODEL integration constants
c IPRT: Unit number for text out put (normally 6 for screen)
c LD: Old single precision flag. No longer used.
c
COMMON/ACCU/EPS,PVSTP,EVSTP,IBG1,IIIBG1,IBG2,IIIBG2,IPHGB,
* NDOCUP(5),STPSZ,EPSMOD
c
c The only elements used/set in mjsanl are:
c IBG1: Flag to turn on the PMODEL detail in CPINT. 6=Yes, 0=No
c EPSMOD: A value of 0.0 turns on all feedthrough.
c
c Set up common block with RMARK values for E1, E2, and ion mode
c
common/lim/rmarkp(4,3)
c
c RMARK(i,j): The currents in fA at particular digital current levels.
c where i: 1=L; 2=M; 3=E1; 4=E2 modes
c and j: 1=digital current at channel 0; 2=128; 3=255.
c
c Time of spectra variables containing mSec from 1977
c
REAL*8 TBEG,TSPEC
LOGICAL LFIRST, REW
c
c Currents in the cups
c
REAL BARR(129,5),BPLT(129,5)
c
c Cup temperatures in Centigrades and current levels in KNTCUR
c
DIMENSION TEMP(3),RMARK(3)
c
c Status word masking variables
c
integer*2 i2m(2)
equivalence (kmask,i2m(1))
c
c Flag for first time around
c
DATA LFIRST/.TRUE./
c
c Status word mask corresponding to KMASK=000F0000
c
data i2m /15,0/
c
c Possible input data tape/file types
c
data idtype/'EDR ','SUM ','SPL ','VGRT'/
c
c Astronomical unit in km
c
DATA AU/1.495979E8/
c
c Clear initially all the currents
c
data BARR/645*0.0/,BPLT/645*0.0/
c
c Cup Normal Unit Vectors
c
DATA CUPN/-.2962,-.1710,-.9397,+.2962,-.1710,-.9397,
* +0.0,.3420,-.9397,.7310,.6816,-.0349/
c
c Offset currents for L long (GS5) at Neptune
c
c The offsets averaged in the solar wind is the file:
c /data3/asz/ions/noise.ave
c The first four lines are the averages, the second four are the
c standard deviations. The four lines correspond to the four
c cups in the order A,B,C,D. The numbers in each line are the
c L mode channels in the order from the lower to higher.
c
c The values used below are selected for offset currents greater than
c standard deviation currents. THESE SHOULD BE REGARDED AS PRELIMINARY!!!
c 3/23/90 rlm
c
data curoff/0.,0.,0.,0.,0.,0.,0.,0., 0., 0.,
+62.,249.,334.,967.,1105.,1349.,
c End A cup
+0.,0.,0.,0.,0.,0.,0.,0.,
+0.,0.,0.,0.,0.,0.,0.,0.,
c End B cup
+0.,0.,0.,0.,0.,0.,0.,0., 0., 0.,
+95.,228.,201.,217.,276.,298.,
c End C cup
+0.,0.,0.,0.,0.,0.,0.,0.,
+0.,0.,0.,0.,0.,0.,0.,0./
c End D cup
c
c==========================================================
c B E G I N E X E C U T I O N
c==========================================================
c
c Initialize trajectory file flag to no traj data available
c
isedr=0
c
c Initialize SC velocity in SC coordinates to 0.0
c
do 1201 i=1,3
sblk(39+i)=0.0
1201 continue
c
c Initialize rotation matrix from SC to
c cylindrical coordinates to the identity matrix
c
do 21 i=1,9
sblk(43+i)=0.0
21 continue
sblk(44)=1.0
sblk(48)=1.0
sblk(52)=1.0
c
c End matrix initialization
c
c Determine the Version of Program
c
c The RUN_TYPE environment variable is set in the RUN command
c executed from the script file. It value is a 3 character string
c which is 'vgr' for cruise mode (mission tapes) and 'spl' for
c encounter mode (spectral tapes). Read this variable to
c determine what is the appropriate mode of operation.
c
call getenv('RUN_TYPE',run_type)
if (run_type .eq. 'vgr') then
isel=0
write(6,*)
write(6,*)'========================================'
write(6,*)' Cruise Mode (Mission tape)'
write(6,*)'========================================'
write(6,*)
else if (run_type .eq. 'spl') then
isel=1
write(6,*)
write(6,*)'========================================'
write(6,*)' Encounter Mode (Spectral tape)'
write(6,*)'========================================'
write(6,*)
else
write(6,*)
write(6,*)'ERROR! Unknown mode of operation.'
write(6,*)' The environment RUN_TYPE is incorrectly set in RUN'
write(6,*)
stop
endif
89 if (isel .eq. 1) go to 201
c
c Select trajectory option in Cruise Mode
c
write (6,*)'Trajectory file selection:'
write (6,*)'0 ==> Do not mount trajectory file'
write (6,*)'1 ==> Use short cruise file (via PHASE3)'
read (5,*,end=89,err=89) isedr
c
c The flag variable ISEDR should have the value 0 for no trajectory
c or 3 for the use of PHASE3.
c
isedr = isedr * 3
go to 202
c
c Select trajectory option in Encounter Mode
c
c
201 continue
write(6,*)'Trajectory file selection:'
write(6,*)'0 ==> Do not mount trajectory file'
write(6,*)'1 ==> Use short encounter file (via PHASE2)'
read(5,*,end=89,err=89) isedr
c
c The flag variable ISEDR should have the value 0 for no trajectory
c or 2 for the use of PHASE2.
c
isedr = isedr * 2
c
c First call to initialize and start reading file
c
202 if (isedr .eq. 2) CALL PHASE2(JTB,IEOD)
if (isedr .eq. 3) CALL PHASE3(JTB,IEOD)
if (isedr .gt. 3) go to 89
if (isedr .lt. 0) go to 89
c
c Setup graphics window
c
call plot( 0., 1., 3)
call plot( 0., .1, 2)
call nupage2 (W,H)
c
c Zero out the common block /BMODEL/
c
ibmodl=0
idsrnn=0
do 7776 i=1,10
ik(i)=0
7776 continue
ik(1)=1
ik(2)=1
ik(3)=1
ik(4)=1
c
c Initialize the common block /CNTRL/
c
icur=0
i5 = 5
iplt = 0
irec = 17
call initio (1,irec,itemp)
c
c Define the plot size in inches
c
W=11.0
H=8.5
c
c Setup runtime info
c
CALL WRTIME(0)
CALL TIMER
CALL HOTIME
c
c Initialize the text output unit number.
c WARNING! This flag flag is used inconsistently. Changing
c its value will redirect only some of the output.
c
IPRT=6
c
c Decimal time initialization
c
TSPEC=0.0
c
c Termination flag: 0=keep current spectra, 1=read next spectra
c 2=terminate, 3=same data
c
iread=1
c
c Differentiate on type of sedr data
c
if (isedr .eq. 0) go to 5998
if (isedr .eq. 2) go to 2
c
c Report mission and sedr tapes connected for cruise mode
c and the beginning and times of these files.
c
call getenv('ft32',ftdata)
call getenv('ft33',ftsedr)
call getenv('SSEDR_TIME',sedrtime)
call getenv('VGR_TIME',vgrtime)
write(6,*)
write (6,83) ftsedr
83 format('Cruise SEDR = ',a64)
write (6,84) sedrtime
84 format('SEDR file covers the time range: ',a24)
write(6,*)
write (6,85) ftdata
85 format('Cruise DATA = ',a64)
write (6,86) vgrtime
86 format('Mission file covers the time range: ',a24)
write(6,*)
goto 111
c
c Warn that no trajectory data is connected
c
5998 continue
write (6,*)'--- No Trajectory Info Available ----'
go to 111
c
c Report encounter tape start and stop times
c
2 CONTINUE
GO TO (5,10,15,20),NUMFOR
20 WRITE(6,6)
6 FORMAT(22x,'*** NEPTUNE DATA DISPLAY ***')
WRITE(6,6000)
6000 FORMAT(' SEDR DATA FILE COVERS: ',5x,
+'SPECTRAL DATA FILE STARTS AT: ')
WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
6001 FORMAT('V2 - ',2I4,I3,I2,' TO '2I4,I3,I2,5x,'V2 - '2I4,I3,I2)
6003 FORMAT('V1 - ',2I4,I3,I2,' TO '2I4,I3,I2,5x,'V1 - '2I4,I3,I2)
WRITE(6,*)
GO TO 111
15 WRITE(6,7)
7 FORMAT(22x,'*** URANUS DATA DISPLAY ***')
WRITE(6,6000)
WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
WRITE(6,*)
GO TO 111
10 WRITE (6,4)
4 FORMAT (22x,'*** SATURN DATA DISPLAY ***')
WRITE(6,6000)
IF (ISCID .EQ. 1) THEN
WRITE(6,6003)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
ELSE
WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
ENDIF
WRITE(6,*)
GO TO 111
5 WRITE(6,8)
8 FORMAT(22x,'*** JUPITER DATA DISPLAY ***')
WRITE(6,6000)
IF (ISCID .EQ.1) THEN
WRITE(6,6003)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
ELSE
WRITE(6,6001)(IST(I),I=1,4),(ISP(I),I=1,4),(JTB(I),I=1,4)
ENDIF
WRITE(6,*)
111 continue
if (iread .ne. 0) go to 1333
write (6,1334) TIMESR
1334 format('Time of last spectrum is ',5a4)
1333 continue
c
c Option to do joint analysis of E1, L, E2 or E2, M, E1
c data sets with PMODEL and all feedthrough code turned on
c
ifull=0
c
c Modified to work on GS5 data with L modes;
c only long modes called - last spectrum read must
c be L-long or M.
c
write (6,3550)
3550 format ('Do joint ion/electron analysis?',
+' 2=M, E1, E2; 1=L, E1, E2; 0=NO')
read (5,*,end=111,err=111) ifull
if (ifull .eq. 0) go to 3551
if ((ifull .eq. 1) .or. (ifull .eq. 2)) go to 1121
go to 111
c
c Set flags appropriately and skip other inquiries
c
1121 isptyp=1
ipltdc=0
ijtl=0
c
c CLEAR BARR; SET ALL ELEMENTS TO ZERO
c
BARR(1,1) = 0.0
CALL MVC (BARR(1,1),BARR(2,1),2576)
c
c VOY has already been set in phase2,3 and stored in /SCID/
c Load this value in V0 in /PLCONS/
c
V0=VOY
c
c Continue execution
c
if (isel .ne. 0) go to 597
go to 777
c
c Select type of data to be examined
c
3551 continue
c
c CLEAR BARR; SET ALL ELEMENTS TO ZERO
c
BARR(1,1) = 0.0
CALL MVC (BARR(1,1),BARR(2,1),2576)
V0=VOY
write (6,3196)
3196 format (1h ,'Select Data Mode:')
write (6,3195)
3195 format (1h ,'1 = Plasma, 2 = DC Return,'
+,' 3 = MODCAL and CURCAL, 4 = All')
read (5,*,END=111,Err=111) isptyp
ijtl=0
c
c Option to plot summed plasma current in DC return
c format; skip analysis
c
6711 continue
WRITE(IPRT,770)
770 FORMAT(1H ,'Integrate Plasma Data? (Analysis',
+' Skipped)')
write (iprt, 769)
769 format(1h ,'2=YES for DC return spectra')
write (iprt, 768)
768 format(1h ,'1=YES for plasma spectra')
write (iprt, 767)
767 format(1h ,'0=NO')
ipltdc = 0
READ(5,*,END=6711,Err=6711) ipltdc
if ((isptyp .ne. 1) .and. (isptyp .ne. 4)) go to 777
c
c Select spectral type if plasma data
c
597 continue
write (6,3199)
3199 format (1h ,'Select Spectral Type:')
write (6,3198)
3198 format (1h ,' 0 = all, 1 = ions only,'
+,' 2 = electrons only')
write (6,3197)
3197 format (1h ,' 3 = L only, 4 = M only, 5 = E1 only,'
+,' 6 = E2 only, 7 = LS, E1S, E2S only')
read (5,*,END=111,Err=111) ijtl
c
c Select time of first data to be viewed
c
777 continue
WRITE (6,100)
100 FORMAT(1H ,'INPUT START TIME: IYR, IDAY, IHR, IMIN, ISEC, IMSEC')
itb(6)=-1
READ (5,*,END=111,Err=111)ITB
if (itb(6) .eq. -1) go to 111
c
c Calculate the number of msec since the beginning of 1977
c to the time specified by the array ITB
c WARNING: This routine does NOT correct for leap years.
c
CALL xTYME(ITB,TBEG)
c
c IF THE DESIRED SPECTRUM IS EARLIER IN TIME THAN THE
c CURRENT ONE, A CALL TO VGRREW IS MADE. THIS WILL REWIND THE TAPE.
c
c phase2 and phase3 takes care of rewinding the traj tapes
c
c First time around TSPEC is zero
c
if (TBEG .ge. TSPEC) go to 774
c
c The routines VGRREW are located in the nxtvgr.f and nxtspl.f
c modules. They have a slightly different format.
c
IF((TBEG.LT.TSPEC) .and. (ntype .ne. 4)) CALL VGRREW
if (ntype .eq. 4) then
REW = .true.
CALL VGRREW(REW)
endif
774 continue
c
c Set read flag to read next spectra
c
IREAD=1
c
c For combination modes, do not allow averaging
c
if ((ijtl .lt. 3) .or. (ijtl .gt. 6)) go to 11
c
c Option for averaging over NS spectra before doing analysis
c (Similar, but not the same, as implementation via CAVE in KNTCUR -
c not implemented in this code). rlm 3/7/90
c
write (6,*) "Input number of spectra to average (0= no ave.)"
read (5,*) ns
nsc=0
do 5003 k=1,512
curave(k)=0.0
ncount(k)=0
5003 continue
GO TO 11
1 CONTINUE
199 continue
WRITE (6,200)
200 FORMAT(1H ,'READ NEXT SPECTRUM? 1=YES, 0=NO, 2=TERMINATE',
+' 3=SAME DATA'/)
write (6,2200)
2200 format (1h ,'"0=NO" allows reset of spectral type selection')
iread = -1
READ(5,*,END=199,Err=199) IREAD
if (iread .eq. -1) go to 199
IF (IREAD .EQ. 1) GO TO 11
IF (IREAD .EQ. 0) GO TO 2
IF (IREAD .EQ. 2) GO TO 1000
IF (IREAD .EQ. 3) GO TO 1101
GO TO 199
11 CONTINUE
CALL NXTMOD (ANS,nJTB,JDATA,JNE,LSTAT,JTLMOD,JCLK,TEMP,IEOD,IL,
*ITYPE)
c
c Check for short versus long data (GS-5)
c
kstat=lstat
call SETJTL( JTLMO, JTLOFF, KSTAT, nJTB, JTLMOD, JCLK)
c
c Returned:
c JTLMO: =JTLMOD unless short spectra then =JTLMOD+4
c JTLOFF: 4 for short, 0 for long modes
c KSTAT: only the low 2 bytes are kept, rest zeroed out.
c
c Copy over time as read by NXTMOD
c
jtb(1)=njtb(1)
jtb(2)=njtb(2)
jtb(3)=njtb(3)
jtb(4)=njtb(4)
jtb(5)=njtb(5)
jtb(6)=njtb(6)
c
c CHECK FOR END OF FILE OR ERROR
c
IF (IEOD .EQ. 0) GO TO 1001
WRITE (6,198)
198 FORMAT ('END OF DATA OR ERROR ENCOUNTERED IN SPECTRAL FILE')
WRITE (6,197)
197 FORMAT ('DO YOU WISH TO CONTINUE WITH CURRENT S/C AND PLANET?')
195 continue
WRITE (6,196)
196 FORMAT ('TYPE 1 FOR YES OR 0 TO TERMINATE')
iend = -1
READ (5,*,END=195,Err=195) IEND
if (iend .eq. -1) go to 195
IF (IEND .EQ. 1) GO TO 111
IF (IEND .EQ. 0) GO TO 1000
GO TO 195
c
c INITIALIZE
c
1001 IF ( .NOT. LFIRST) GO TO 1101
LFIRST=.FALSE.
CALL RUNBEG
CALL PLSBEG
CALL CONRD
c
c Initialize values in the common block /SCALE/
c
YB=.1
IDEC=6
XB=0.
c
c Put data tape/file type into /TPTYPE/
c
jtype=itype
c
c Find corresponding flag number and put it into /TPTYPE/
c NTYPE=1 - EDR; 2 - SUM; 3 - SPL; 4 - VGRT
c
do 505 ii=1,4
if (itype .eq. idtype(ii)) nntype=ii
505 continue
ntype=nntype
c
c End initialization
c
1101 continue
c
c Reset JTLMO, JTLOFF and KSTAT
c
call GETJTL( JTLMO, JTLOFF, KSTAT)
c
c Check type for full up analysis and reject unselected modes
c
if ((ifull .eq. 1) .and. (jtlmo .eq. 2)) go to 11
if ((ifull .eq. 2) .and. (jtlmo .eq. 1)) go to 11
c
c DECODE STATUS TO DETERMINE INSTRUMENT MODE
c
c Keep only one byte of the status word
c
IPLS=IAND(LSTAT,KMASK)
c
c Move that one byte to the lowest position
c
IPLS=IPLS/2**16
c
c The meaning of the lowest byte of the status word:
c Status byte = XX Y Z
c 84 2 1
c
c Z = 1/0 = Modulator On/Off
c Y = 1/0 = Plasma / DC Return mode
c XX = 00 = PLS reverse grid
c 01 = Modulator voltage monitor
c 10 = Current calibration
c 11 = PLS normal grid
c
c In case all modes were selected, do not check the status byte
c
if (isptyp .eq. 4) go to 900
c
c IDC - MODULATOR ON (Normal .OR. Reversed)
c
IF ((IPLS .EQ. 13 .OR. IPLS .EQ. 1) .and. isptyp .eq. 2) GO TO 112
c
c IDC - MODULATOR OFF (Normal .OR. Reversed)
c
IF ((IPLS .EQ. 12 .OR. IPLS .EQ. 0) .and. isptyp .eq. 2) go to 112
c
c PLASMA DATA - MODULATOR ON (Normal .OR. Reversed)
c
IF ((IPLS .EQ. 15 .OR. IPLS .EQ. 3) .and. isptyp .eq. 1) GO TO 900
c
c PLASMA DATA - MODULATOR OFF (Normal .OR. Reversed)
c
IF ((IPLS .EQ. 14 .OR. IPLS .EQ. 2) .and. isptyp .eq. 1) GO TO 900
c
c Modulator Voltage Monitor (MVM / MODCAL)
c (Voltage Monitor On + Modulator On)
c
IF (IPLS .EQ. 7 .and. isptyp .eq. 3) go to 112
c
c Current Calibration (CURCAL)
c Calibrate Measurement Chain (Modulator Off) = A
c
IF (IPLS .EQ. 10 .and. isptyp .eq. 3) go to 112
c
c Background Correction (Modulator On) = B
c
IF (IPLS .EQ. 11 .and. isptyp .eq. 3) go to 112
c
c Illegal Status: IPLS = 4, 5, 6, 8, 9
c
GO TO 11
900 CONTINUE
c
c CONTINUE ONLY for selected modes if plasma data
c
if (jtlmo .gt. 4) go to 3193
if (ijtl .eq. 0) go to 112
idiff=ijtl*1.-jtlmod*1.
if (idiff .eq. 2) go to 112
if (ijtl .eq. 1 .and. jtlmod .lt. 3) go to 112
if (ijtl .eq. 2 .and. ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)))
+go to 112
go to 11
c
c Only accept short modes if it was specificaly requested
c
3193 continue
if ((ijtl .eq. 7) .and. (jtlmo .ge. 5)) go to 112
go to 11
c
c SEARCH FOR FIRST SPECTRUM AFTER START TIME
c
112 continue
jtlold=jtlmod
mp=jtlmod-1
me=4-jtlmod
c
c Store current spectra time in /RTIM/
c
DO 75 N=1,6
RCUR(N)=JTB(N)
KTB(N)=JTB(N)
75 CONTINUE
c
c Convert current time to msec since 1977
c
CALL xTYME(KTB,TSPEC)
c
c If current time is before the requested time, read next spectra
c
IF (TSPEC .LT. TBEG) GO TO 11
c
c DATA NUMBER TO ENGINEERING UNIT (FEMTOAMPERE) CONVERSION (NO FILTER)
c (KNTCUR does no decoding for MVM; uses stored lookup table for
c DC returns - apparently from Anton's calibration data for -10 de C;
c PLASMA and CURCAL currents converted based on current status word
c for amplifier gain and capacitor and on JCLK for current integration
c time). Negative JNE turns off vgranl filtering in KNTCUR.
c
CALL KNTCUR(ANS,JTB,JDATA,-JNE,LSTAT,JTLMOD,JCLK,TEMP,CURRNT,KSTAT
*,IXCEL,CNOISE,RMARK)
c
c INPUT: The standard calling sequence described in NXTMOD
c OUTPUT:
c CURRNT(512): The measured currents in each channel (fA)
c if -1 - suspect data; -2 - data not received
c KSTAT: Best guess status word
c IXCEL: Quality flag: 1 = Ok. If larger positive, it is the
c number of suspect channels. If negative, some channels
c are saturated.
c CNOISE: The noise level for the current instrument settings
c RMARK(3): The current levels for digital currents 0, 128, 255.
c
c Skip averaging if it was not asked for
c
if (ns .eq. 0) go to 5005
c
c Collect spectra for averaging. Accept only currents above
c the minimum threshold.
c
nsc=nsc+1
do 5000 k=1,4*jne
c if (currnt(k) .gt. rmark(1)) ncount(k)=ncount(k)+1
if (currnt(k) .ge. rmark(1)) ncount(k)=ncount(k)+1
c if (currnt(k) .gt. rmark(1)) curave(k)=curave(k)+currnt(k)
if (currnt(k) .ge. rmark(1)) curave(k)=curave(k)+currnt(k)
5000 continue
jtlold=jtlmod
if (nsc .lt. ns) go to 11
c
c Once collected the right number of spectra, divide by count
c
do 5001 k=1,4*jne
if (ncount(k) .ne. 0) currnt(k)=curave(k)/(ncount(k)*1.)
5001 continue
c
c Reset variables used for averaging
c
nsc=0
do 5002 k=1,4*jne
ncount(k)=0
curave(k)=0.0
5002 continue
5005 continue
c
c Transfer threshold information to common block /LIM/
c
do 700 j=1,3
rmarkp(jtlmod,j)=rmark(j)
700 continue
c
c For Cruise data and SUM or long SPL tapes, load
c initial values from answer array; key on
c ntype = 4 (passed through common block TPTYPE)
c
c Check if Summary tape
c
if (nntype .eq. 2) ntype=4
c
c Check if long spectral tape
c
if ((nntype.eq.3) .and. (ans(38).lt.1000.0) .and.
& (ans(38) .ne. 0.0)) ntype=4
775 continue
c
c Get trajectory information for the requested time
c
if (isedr .eq. 3) CALL PHASE3(JTB,IEOD)
if (isedr .eq. 2) CALL PHASE2(JTB,IEOD)
c
c Do not clear barr each time at this point
c
c TRANSFER CURRENTS TO BARR
c Put E1 or E2 currents into D-cup
c
if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 1450
c
c Place the # of channels in a cup into /FULL/
c
jnei=jne
DO 495 I=1,4
DO 595 K=1,JNE
BARR(K,I)=CURRNT(K+JNE*(I-1))
c
c Also put ion currents into CURNEW in /FULL/
c
curnew(K+JNE*(I-1))=CURRNT(K+JNE*(I-1))
595 CONTINUE
495 CONTINUE
go to 1350
c
c Move electron currents into D cup slot in BARR if IFULL=0
c
1450 continue
if ((ifull .eq. 1) .or. (ifull .eq. 2)) go to 3553
DO 3595 K=1,JNE
BARR(K,4)=CURRNT(K)
BPLT(K,4)=BARR(K,4)
3595 CONTINUE
go to 3554
c
c Put E1 electrons into Barr(1-16,5) for full up analysis
c Put E2 electrons into Barr(17-32,5)
c
3553 continue
do 3555 k=1,16
j=(jtlmod-3)*16+k
barr(j,5)=currnt(k)
3555 continue
c
c Move electrons into D cup slot in CURRNT
c for both full and not full up analysis
c
3554 continue
DO 3495 k=1,JNE
CURRNT(K+JNE*3)=BARR(K,4)
j=(jtlmod-3)*16+k
curnew(j+jnei*4)=barr(j,5)
3495 continue
1350 continue
c
c Place the # of channels per cup into /SCID/ and /SCALE/
c
IME=JNE
XSCALE=JNE
c
c WRITE OUT TIME OF CURRENT SPECTRUM STORED in KTB
c
CALL YDOYHM(KTB(1),KTB(2),KTB(3),KTB(4),KTB(5)+KTB(6)/1000.)
c
c YDOYHM fills the /CDATES/ common block including the 20 character
c string TIMESR containing the current spectrum time.
c
WRITE(IPRT,38) TIMESR
38 FORMAT('TIME OF CURRENT SPECTRUM IS ',5A4)
IF (JNE .EQ. 16 .AND. JTLMOD .EQ. 1) WRITE (IPRT,747)
747 FORMAT('CURRENT SPECTRUM IS AN L-MODE')
IF (JNE .EQ. 16 .AND. JTLMOD .EQ. 4) WRITE (IPRT,47)
47 FORMAT('CURRENT SPECTRUM IS AN E2-MODE')
IF (JNE .EQ. 16 .AND. JTLMOD .EQ. 3) WRITE (IPRT,746)
746 FORMAT('CURRENT SPECTRUM IS AN E1-MODE')
IF (JNE .EQ. 128) WRITE (IPRT,48)
48 FORMAT('CURRENT SPECTRUM IS AN M-MODE')
c
c For Titan Proposal, if not at Saturn, skip this section
c
if (numfor .ne. 2) go to 8889
c
c Assume Titan radius = 2575 km and calculate the
c distance of the S/C from the cloud tops.
c
titalt=sblk(55)-2575
write (6,8888) sblk(55),titalt
8888 format (' Distance to Titan center =',f9.2,
+' = altitude of ',f9.2)
c
8889 continue
c
c Find B field projections along cup normals
c
if (ntype .eq. 4) WRITE(6,5666)
5666 FORMAT('PROJECTION OF Magnetic Field ALONG CUP NORMALS',
* ' (Gamma)')
if (ntype .eq. 4) WRITE(6,5667)
5667 FORMAT(1H ,5X,5X,'BA',9X,'BB',9X,'BC',9X,'BD',9x,'B mag')
DO 5671 I=1,4
BCUPN(I)=-(ans(5)*CUPN(1,I)+ans(6)*CUPN(2,I)+
* ans(7)*CUPN(3,I))
5671 CONTINUE
if (ntype .eq. 4) WRITE(6,5672)(BCUPN(I),I=1,4),ans(8)
5672 FORMAT(3X,5(2x,F9.3)/)
c
c End insert for B field projection
c
c WRITE OUT TYPE OF DATA
c
IF (IPLS .EQ. 3) WRITE (6,910)
910 FORMAT ('MODULATOR ON, PLASMA, ** REVERSED SUPPRESSOR **')
IF (IPLS .EQ. 15) WRITE (6,911)
911 FORMAT ('MODULATOR ON, PLASMA, NORMAL SUPPRESSOR')
IF (IPLS .EQ. 2) WRITE (6,3910)
3910 FORMAT ('MODULATOR *OFF*, PLASMA, ** REVERSED SUPPRESSOR **')
IF (IPLS .EQ. 14) WRITE (6,3911)
3911 FORMAT ('MODULATOR *OFF*, PLASMA, NORMAL SUPPRESSOR')
IF (IPLS .EQ. 13) WRITE (6,912)
912 FORMAT ('MODULATOR ON, DC RET, NORMAL SUPPRESSOR')
IF (IPLS .EQ. 1) WRITE (6,913)
913 FORMAT ('MODULATOR ON, DC RET, ** REVERSED SUPPRESSOR **')
IF (IPLS .EQ. 12) WRITE (6,914)
914 FORMAT ('MODULATOR *OFF*, DC RET, NORMAL SUPPRESSOR')
IF (IPLS .EQ. 0) WRITE (6,915)
915 FORMAT ('MODULATOR *OFF*, DC RET, ** REVERSED SUPPRESSOR **')
If (ipls .eq. 7) write (6,916)
916 format ('Modulator Voltage Monitor Test')
if (ipls .eq. 10) write (6,917)
917 format ('Current Calibration Test - Modulator Off')
if (ipls .eq. 11) write (6,918)
918 format ('Current Calibration - Mod. On - Bckgnd. Cal.')
101 CONTINUE
c
c Voltage or current calibration
c
if ((ipls .eq. 7) .or. (ipls .eq. 10) .or.
+(ipls .eq. 11)) go to 6695
c
c Anything but plasma mode with normal or reverse grid with
c either the modulator on or off.
c
if ((ipls .ne. 3) .and. (ipls .ne. 15) .and.
+(ipls .ne. 2) .and. (ipls .ne. 14)) go to 995
c
c For electron modes, skip the next max response section
c
if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 507
c ----------------------------------------------------------------------
c
c FIND CUP WITH MAXIMUM RESPONSE TO COROTATING, COLD BEAM OF PROTONS
c
AMAXH=-1.
IPH=0
DO 500 K=1,4
IF (AMAXH .GE. SBLK(8+K)) GO TO 500
AMAXH=SBLK(8+K)
IPH=K
500 CONTINUE
ANS(20)=0.
IF (AMAXH .LE. 0. .OR. SBLK(IPH+12) .LE. 0.) GO TO 507
SUM=0.
DO 501 I=1,JNE
SUM=SUM+BARR(I,IPH)
501 CONTINUE
RCUP=9.264E-4
IF (IPH .EQ. 4) RCUP=1.110E-3
ANS(20)=RCUP*SUM/(AMAXH*SBLK(IPH+12))
ANS(21)=IPH
c
c Get density from mission tape if it is not the filll value of `preset`
c
if ((ntype .eq. 4) .and. (ans(38) .lt. 1.e6)) ans(20)=ans(38)
c --------------------------------------------------------------------
507 CONTINUE
c
c ASK IF TO PLOT REAL CURRENTS
c
670 continue
WRITE(IPRT,70)
70 FORMAT('PLOT DATA? 1=YES, 0=NO')
iplt = -1
READ(5,*,END=670,Err=670) IPLT
if (iplt .eq. -1) go to 670
IF (IPLT .EQ. 1) GO TO 71
if (iplt .eq. 0 .and. ipltdc .eq. 1) go to 6695
IF (IPLT .EQ. 0) GO TO 60
GO TO 670
c
c Find thresholds and limits
c
71 CONTINUE
thresh=rmark(1)
WRITE (IPRT,53) RMARK(1)
53 FORMAT('THRESHOLD IS ',1PE12.4,' FEMTOAMPS.')
WRITE (IPRT,153) RMARK(3)
153 FORMAT('SATURATION IS ',1PE12.4,' FEMTOAMPS.'/)
WRITE (6,154)
154 FORMAT ('MINIMUM AND MAXIMUM CURRENTS (IN FEMTOAMPS) ARE:')
if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 3154
DO 95 J=1,4
IS=JNE*(J-1)+1
CALL LIMITS (CURRNT(IS),IME,AMIN,AMAX,IMIN,IMAX)
WRITE (IPRT,50) AMIN,AMAX,J
50 FORMAT(2(1X,1PE12.3),' IN CUP ',I5)
95 continue
go to 995
3154 continue
CALL LIMITS (BPLT(1,4),IME,AMIN,AMAX,IMIN,IMAX)
WRITE (IPRT,350) AMIN,AMAX
350 FORMAT(2(1X,1PE12.3))
995 continue
write (iprt,*)
WRITE (IPRT,430)
430 FORMAT('WRITE OUT ALL CURRENTS? 1=YES, 0=NO')
irite = -1
READ (5,*,END=995,Err=995) IRITE
if (irite .eq. -1) go to 995
IF (IRITE .le. 0) GO TO 6697
IF (IRITE .EQ. 1) GO TO 1695
GO TO 995
1695 continue
if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 4695
c
c Write Out Plasma Currents for L and M Modes
c
1696 continue
c
c For M mode only give a way out
c
IF (MP .EQ. 1) THEN
WRITE (6,940)
READ (5,*,END=1696,Err=1696) IAN
IF (IAN .EQ. 0) GO TO 6697
ENDIF
DO 1941 J=1,4
WRITE(IPRT,440) J
440 FORMAT(' VOYAGER DATA - CURRENTS (FA) CUP=',I3/)
IM = 16
IF (MP.EQ.1 .AND. J.NE.5) IM=128
WRITE (IPRT,450) (BARR(I,J),I=1,IM)
450 FORMAT(17(1X,1P8E10.3/))
1940 continue
if ( (j .ne. 2) .or. (mp .ne. 1)) go to 1941
WRITE (6,940)
940 FORMAT ('TYPE 1 TO CONTINUE ')
READ (5,*,END=1940,Err=1940) IAN
1941 continue
1942 continue
WRITE (6,940)
READ (5,*,END=1942,Err=1942) IAN
c
c Finished option for writing out currents
c Now do plots; use PLTVGR for normal ions, PLTELE for normal
c electrons; and, PLTFULL for full analysis
c
6697 continue
if (ifull .ne. 0) go to 3556
if ((jtlmod .eq. 3) .or. (jtlmod .eq. 4)) go to 6698
c
c Go to subroutine PLTVGR to actually do plotting
c
if ((ipls .ne. 3) .and. (ipls .ne. 15) .and.
+(ipls .ne. 2) .and. (ipls .ne. 14)) go to 6695
c
c Plot only for plasma mode
c
CALL PLTVGR(BARR,ANS,RMARK)
c
c Case of DC style plots and no analysis for plasma data
c
if (ipltdc .eq. 1) go to 6695
go to 142
c
c Write Out Plasma Currents for E1 and E2 Modes
c
4695 continue
IF (JTLMOD .EQ. 3) WRITE(IPRT,6440)
6440 FORMAT(' VOYAGER DATA - CURRENTS (FA) E1 MODE'/)
IF (JTLMOD .EQ. 4) WRITE(IPRT,4440)
4440 FORMAT(' VOYAGER DATA - CURRENTS (FA) E2 MODE'/)
WRITE (IPRT,6450) (BARR(I,4),I=1,16)
6450 FORMAT(17(1X,1P8E10.3/))
c
c Go to subroutine PLTELE to actually do plotting
c
6698 continue
if ((ipls .ne. 3) .and. (ipls .ne. 15) .and.
+(ipls .ne. 2) .and. (ipls .ne. 14)) go to 6695
c
c Plot only for plasma mode
c
CALL PLTELE(BPLT,ANS,JTLMOD)
c
c Case of DC style plots and no analysis for plasma data
c
if (ipltdc .eq. 1) go to 6695
go to 142
6695 CONTINUE
c
c Plot out and/or analyze calibration, DC, return, or PLASMA MOD OFF
c in subroutine CALANL
c
call calanl(ANS,JTB,JDATA,JNE,LSTAT,JTLMOD,JCLK,TEMP,IEOD,IL,
*ITYPE,ipls,ipltdc,RMARK,thresh)
c
c Return to query "READ NEXT SPECTRUM?"
c
go to 1
c
c OPTION TO DO NONLINEAR LEAST SQUARES ANALYSIS on plasma data
c (data may be from either normal or reversed suppressor modes)
c
c N.B. If DC plots for plasma data are selected (ipltdc=1), then
c the logic is routed through calanl and no analysis is
c performed on this spectrum/spectral set.
c
c
c Plot call for full analysis scheme
c
3556 continue
call PLTFULL(barr, ans)
142 continue
i5 = 5
60 CONTINUE
WRITE (IPRT,42)
42 FORMAT('DO NONLINEAR LEAST SQUARES ANALYSIS WITH CURRENT SPECT',
+'RA? INPUT NUMBER of ITERATIONS.')
write (iprt,1142)
1142 format ('A negative number uses unit 15 as a remote input',
+' file')
if (ipls .eq. 3) write (6,242)
242 format (' WARNING - REVERSED SUPPRESSOR CONFIGURATION')
ifit = -999
read ( i5,*,END=142,Err=142) IFIT
if (ifit .eq. -999) go to 142
nrep = iabs(ifit)
IF (IFIT .EQ. 0) GO TO 1
IF (IFIT .gt. 0) GO TO 8999
i5 = 15
call initioco( 1, i5, itemp)
C
C SELECT PRECISION
C
8999 continue
if (ifull .eq. 0) go to 3558
write (6,3559)
3559 format ('Do you really want to do joint ion/electron analysis '
+'with the current spectra? 1=YES, 0=NO')
read (i5,*, end=8999,err=8999) if2
if (if2 .eq. 0) go to 111
if (if2 .ne. 1) go to 8999
c
c Really do want to continue, so set up all defaults
c for PMODEL integration steps
c
c Include unity response XJC option in as check out
c of code; no cross talk between ion and electron modes
c
3761 continue
write(6,*)'Do you want to use PMODEL with ion/electron ',
& 'interaction (1)'
write(6,*)'or unity response with no ion/electron interaction ',
& '(2)'
read(i5,*,end=3761,err=3761) icplt
if (icplt .eq. 2) icplt=3
if ((icplt .ne. 1) .and. (icplt .ne. 3)) go to 3761
c
c Option for intermediate output
c
6659 write (6,6660)
6660 format ('Turn on PMODEL detail in CPINT? 1=Yes, 0=No')
read (i5,*,end=6659,err=6659) ibg1
if (ibg1 .eq. 1) ibg1=6
if ((ibg1 .ne. 0) .and. (ibg1 .ne. 6)) go to 6659
c
c Set old single precision flag
c
ld=.false.
c
c Set nominal integration parameters for PMODEL
c
c EPSMOD=0.0 turns on all feedthrough
c R1=1.4 sets break from linear to logarithmic integration at 1.4
c thermal widths
c R2=3.0 sets break from logarithmic to inverse integration at 3.0
c thermal widths
c Empirically, these values seem to do the integration rapidly and
c accurately
c
c RSTEP of 7 does adequate job for E2 feedthrough into M
c RSTEP > 4 is probably necessary for a reasonable job
c
c RSTEP=7.0
c rstep=3.5
c
c Species into nominal mode
c
c The following values were used for fit shown at 3rd Brice
c meeting in Lindau, FRG, Oct. 10-13, 1988. Fit was to
c two O+ and three e- Maxwellians. Both n and w were fit for
c each (but NOT the velocity vector). Four iteractions gave
c convergence; 16 cpu minutes per iteration on PLASMA using the
c floating point accelerator board.
c
cc rstl=5.0
cc rstm=0.1
cc rste1=4.0
cc rste2=0.7
cc rstle=2.0
cc rstme=0.5
cc rste1i=0.2
cc rste2i=0.08
c
c New defaults (test with V1 Saturn data 10/26/88)
c
rstl=6.0
c
c M mode parameters neither checked nor verified
c 10/26/88 rlm
c
rstm=0.5
rste1=1.0
rste2=0.5
rstle=2.0
rstme=0.5
rste1i=0.5
rste2i=0.1
3662 continue
write (6,3661) rstl,rstm,rste1,rste2
3661 format ('Default for RSTEPs (direct):',4(1x,f6.3),
+' for L, M,',' E1, E2, respectively.')
3665 write (6,3664)
3664 format ('Do you wish to override? 0=NO')
read (i5,*,end=3662, err=3662) ior
if (ior .eq. 0) go to 3672
c
c Change defaults
c
3666 continue
write (6,3668)
3668 format ('Input new values (direct)')
read (i5,*,end=3666,err=3666) rstl,rstm,rste1,rste2
write (6,3663) rstl,rstm,rste1,rste2
3663 format ('RSTEPs = ',4(1x,f6.3))
c
c Species into feedthrough mode
c
3672 continue
write (6,3671) rstle,rstme,rste1i,rste2i
3671 format ('Default for RSTEPs (fdthr):',4(1x,f6.3),
+' for L, M,',' E1, E2, respectively.')
3675 write (6,3674)
3674 format ('Do you wish to override? 0=NO')
read (i5,*,end=3672, err=3672) ior
if (ior .eq. 0) go to 3682
c
c Change defaults
c
3676 continue
write (6,3678)
3678 format ('Input new values (fdthr)')
read (i5,*,end=3676,err=3676) rstle,rstme,rste1i,rste2i
write (6,3673) rstle,rstme,rste1i,rste2i
3673 format ('RSTEPs = ',4(1x,f6.3))
c
c End RSTEP specification
c
3682 continue
R1=0.7
R2=3.0
EPSMOD=0.0
3045 continue
WRITE (6,3040) R1,R2
3040 FORMAT(' CURRENT PMODEL INTEGRATION CONSTANT SETTINGS ARE:'/,
1' R1 = ',F5.2,' R2 = ',F5.2/,
2' DO YOU WISH TO OVERRIDE? 1=YES, 0=NO')
iover = -1
read ( i5,*,END=3045,Err=3045) IOVER
if (iover .eq. -1) go to 3045
IF (IOVER .EQ. 0) GO TO 3020
IF (IOVER .EQ. 1) GO TO 3142
GO TO 3045
3142 continue
WRITE (6,3042)
3042 FORMAT('R1 = ')
read ( i5,*,END=3142,Err=3142) R1
3143 continue
WRITE (6,3043)
3043 FORMAT('R2 = ')
read ( i5,*,END=3143,Err=3143) R2
GO TO 3045
3020 continue
write(6,3560)
3560 format ('Start joint ion/electron analysis sequence')
c
c Call XFULL to start analysis
c
go to 4021
c
c Normal analysis/options package follows
c
3558 if (jtlmod .lt. 3) go to 8998
write (iprt,243)
243 format ('Analysis code not implemented for E1 or E2')
go to 1
8998 continue
c
c No differentiation between single and double precision
c on SUNs - this is a leftover from O.I.T./GSFC IBM
c 3/10/89 - rlm
c
ld=.false.
c
c SELECT SIMULATION OPTION
c
go to 8996
8997 continue
i5 = 5
8996 continue
WRITE (IPRT,41)
41 FORMAT (" SELECT FITTING OPTION:",
| " if negative calculate current " )
WRITE (6,43)
43 FORMAT (8X,'0 ==> DO NOT DO FIT')
WRITE (6,44)
44 FORMAT (8X,'1 ==> FIT USING PMODEL')
c WRITE (6,45)
c 45 FORMAT (1H ,8X,'2 ==> FIT USING XJC')
WRITE (6,52)
52 FORMAT (8X,'3 ==> FIT USING XJC (UNITY RESPONCE)')
WRITE (6,46)
46 FORMAT (8X,'4 ==> FIT USING ABSO (L OR M MODES - LABCUR,LDCUR)')
WRITE (6,55)
55 FORMAT (8X,'5 ==> FIT USING ABSOM (M-MODE ONLY - ABCCUR,DCUR)')
WRITE (6,123)
123 FORMAT (8X,'6 ==> FIT USING ABSO (BIMAXWELLIAN VERSION)')
WRITE (6,51)
51 FORMAT (8X,'7 ==> FIT USING COLD BEAM (CUPINT,DCPINT)')
c WRITE (6,49)
c 49 FORMAT (8X,'8 ==> FIT USING RECTANGULAR APPROX (XCUR)')
WRITE (6,149)
149 FORMAT (8X,'9 ==> FIT USING VERY COLD BEAM (CUPEXP)')
c WRITE (6,163)
c 163 FORMAT (1H ,8X,'10 ==> FIT USING FAST VERSION OF XCUR (XTEST)')
WRITE (6,263)
263 FORMAT (5X,'GT 10 ==> TERMINATE EXECUTION ',$)
icplt = -999
read ( i5,*,END=8997,Err=8997) ICPLT
if (icplt .eq. -999) go to 8997
IF (ICPLT .EQ. 0) GO TO 1
icur = icplt
icplt = iabs(icplt)
IF (ICPLT .GT. 10) GO TO 1000
c
c CONSISTENCY CHECK
c
if ( (icplt .eq. 2) .or. (icplt .eq. 8) .or. (icplt .eq. 10))
| go to 8997
IF (ICPLT .GE. 7) GO TO 4020
IF (ICPLT .EQ. 1 .AND. JTLMOD .EQ. 2) GO TO 4000
IF (ICPLT .EQ. 5 .AND. JTLMOD .EQ. 1) GO TO 4010
IF (ICPLT .EQ. 1) GO TO 4030
GO TO 4020
4000 continue
WRITE (6,4001)
4001 FORMAT (1H /,' ** WARNING ** -- "PMODEL" SELECTED FOR M-MODE ANAL
1YSIS. LARGE AMOUNTS OF CPU'/' TIME MAY BE USED. DO YOU WISH TO CON
2TINUE? 1=YES, 0=NO'/)
icont = -1
read ( i5,*,END=4000,Err=4000) ICONT
if ( icont .eq. -1) go to 4000
IF (ICONT .EQ. 0) GO TO 60
IF (ICONT .EQ. 1) GO TO 4030
GO TO 4000
4010 WRITE (6,4011)
4011 FORMAT (1H /,' "ABSOM" CAN SIMULATE M-MODE DATA ONLY'//)
GO TO 60
c
c SET PMODEL INTEGRATION CONSTANTS
c
4030 CONTINUE
RSTEP=7.0
R1=1.4
R2=3.0
EPSMOD=0.0
4045 continue
WRITE (6,4040) RSTEP,R1,R2,EPSMOD
4040 FORMAT(1H ,' CURRENT PMODEL INTEGRATION CONSTANT SETTINGS ARE:'//
1' RSTEP = ',F5.2/
1' R1 = ',F5.2/' R2 = ',F5.2/' EPSMOD = ',F5.2//
2' DO YOU WISH TO OVERRIDE? 1=YES, 0=NO')
iover = -1
read ( i5,*,END=4045,Err=4045) IOVER
if (iover .eq. -1) go to 4045
IF (IOVER .EQ. 0) GO TO 4020
IF (IOVER .EQ. 1) GO TO 4121
GO TO 4045
4121 continue
WRITE (6,4041)
4041 FORMAT('RSTEP = ')
read ( i5,*,END=4121,Err=4121) RSTEP
4142 continue
WRITE (6,4042)
4042 FORMAT('R1 = ')
read ( i5,*,END=4142,Err=4142) R1
4143 continue
WRITE (6,4043)
4043 FORMAT('R2 = ')
read ( i5,*,END=4143,Err=4143) R2
4144 continue
WRITE (6,4044)
4044 FORMAT('EPSMOD = ')
read ( i5,*,END=4144,Err=4144) EPSMOD
GO TO 4045
c
c Do special full up analysis
c
4021 continue
if (jtlmod .gt. 2) go to 61
call xfull(ANS,JNE,JTLMOD,CURRNT,IXCEL,CNOISE,RMARK)
go to 60
61 continue
write (6,62)
62 format (' Last spectrum read must be L or M')
C
C DO NONLINEAR LEAST SQUARES ANALYSIS
C
4020 CONTINUE
icrt=0
if (jtlmod .ne. 1) go to 8112
if (iread .eq. 3) go to 8112
write (6,8111)
8111 format (' Correct L mode currents for noise? 1=Yes, 0=No')
read (i5,*) icrt
if (icrt .ne. 1) go to 8112
do 8113 ikj=1,64
currnt(ikj)=currnt(ikj)-curoff(ikj)
if (currnt(ikj) .le. rmark(1)) currnt(ikj)=rmark(1)
8113 continue
8112 continue
CALL XANAL(ANS,JNE,JTLMOD,CURRNT,IXCEL,CNOISE,RMARK)
GO TO 60
1000 CONTINUE
RETURN
END