* 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