C 06/25/86 15:53:58 06/25/86 15:53:05 x380300.MJS C********************************************************************** C***** THIS ROUTINE WAS MODIFIED BY MRS ON 10/8/85 TO INCREASE ***** C***** THE NUMBER OF POSSIBLE SPECIES FROM 8 TO 15 AND THE NUMBER ***** C***** OF FIT PARAMETERS FROM 18 TO 50. ***** C********************************************************************** c***** modified Jann 20, 1987 by g. s. Gordon Jr. c***** to decrease effect if input typos c********************************************************************** SUBROUTINE XANAL (RES,NE,MD,CURR,IXCEL,CNOISE,RMARK) C C ---------- C | MJSANL | C ---------- C | C ------------------------------------------------------ C | | | | | | C ---------- --------- --------- --------- --------- --------- C | NXTMOD | | KNTCUR | | GETFLD | | PHASE | |PLTVGR | | XANAL | C ---------- --------- --------- --------- --------- --------- C | C | C ---------------------------------------------------- C | | | | | C --------- --------- ---------- ---------- ---------- C | XPARM | | XPACK | | XMJSFT | | XUNPCK | | PLTSIM | <--- C --------- --------- ---------- ---------- ---------- | C | | C |------------------------------------ C --------------- C | | C --------- ---------- C | MJSINV | | XFNCDV | C ---------- ---------- C | C ---------- ( ---------- ) C | XSIMUL |==( | DSIMUL | ) C ---------- ( ---------- ) C | C | C -------------------------------------------------------- C | | | | | | C ---------- ---------- --------- ---------- -------- | C | PMODEL | | LABCUR | | LDCUR | | ABCCUR | | DCUR | | C ---------- ---------- --------- ---------- -------- | C | C | C -------------------------------------------------------- C | | | | | | C ---------- ---------- ------- -------- ---------- | C | CUPINT | | DCPINT | | XJC | | XCUR | | CUPEXP | | C ---------- ---------- ------- -------- ---------- | C | C -------------------------------------------------------- C | | C --------- --------- C | LABBI | | XTEST | C --------- --------- C C RLM 4/17/84 C IMPLICIT INTEGER*2 (I-N) c REAL*8 xxxDT(6) /'xxDATExx','06/25/86','15:53:05', c | ' XANAL',' FORTRAN','MJS '/ COMMON /PLCONS/KDSRN(15),ISC 1 ,VJNE(16,2),DVDE(16,2),CJN(16,2),VOLTSC,J1E(2),J2E(2),FON(2), 1 J3E(2),J4E(2),EN(2),SN(2) 2 ,HAC(3),SIGNOI,DCRTN(256) 3 ,VJNL(16),VJWL(17),VJNM(128),VJWM(129) 4,II,IEND,EPS,FLAM,QK(4,2),NK(4,2),SK(4,2),WK(4,2),VK(3,4), 4 DELTAL(2),IPQF(2,2),IPRF(2,2),IPQ(4,2),IPR(4,2) 5 ,CN(4,3),COSB(3),COS2B(3),SIN2B(3) 6 ,CALCUR(12),CLOCK(3),CAP(2),GAIN(2),V0,ALPHA 7 ,NQ(4), MASS(4),KEY(4),LANIS(4),NCOM(2),NVAR(2),FNOIS 8 ,CSIG,FSIG COMMON/CHIW/WT2(512),CH2 COMMON/CCA/LSYS,LMAG,LPLS C C SET UP TO FIT AS MANY AS 15 DIFFERENT IONIC SPECIES SIMULTANEOUSLY C COMMON/LFIT/LV(3,15),LW(15),LM,LTH,LDEN,LN(15) COMMON/MAG/B(3),ALAM COMMON/CNTRL/IPLT COMMON/NLLS/ICPLT,MP,ME,RSTEP,R1,R2,IPRT,LD COMMON/TIMExx/ITT,IAA INTEGER*4 ITT,IAA,ICPU ,JCPU COMMON/CDATES/SPARA(17),TIMESR(5) COMMON/DAY/DATE,TIME,IQUAL COMMON/ABPARM/MODE COMMON/ROTMAX/TCYSSC(9),TSCRTN(9) DIMENSION VRTN(3,15),VSWRTN(3) DIMENSION SWV(3,15),SWVSV(3),VCUPN(4,15),CUPN(3,4) c DIMENSION LLMM(2) DIMENSION XDEN(15),XVR(15),XVPHI(15),XVZ(15),XWPER(15),XWPAR(15) COMMON/SPEC/SBLK(60) C DATA LLMM/' M',' L'/,SSSS/'SSSS'/,GGGG/'GGGG'/,FFFF/'FFFF'/, C 1EEEE/'EEEE'/ COMMON/MID/NUMFOR LOGICAL LD INTEGER*4 NUMFOR,JSETS(4) C C COMMON BLOCK CONTAINS ALL SIMULATED CURRENTS FROM LAST ITERATION C COMMON/SIMCUR/BTOT(129,5,6,15) C REAL*4 BARR(129,5) LOGICAL*4 LQUAL LOGICAL*4 LV,LW,LM,LTH,LDEN,LN LOGICAL LSYS(32),LMAG(32),LPLS(32),LANIS,LALPHI INTEGER*4 NCMP,MP,ME,ICPLT,MODE,IPRT INTEGER*4 NCOMP,J,N,M INTEGER*4 IPLT DIMENSION RES(150),RMARK(3),VEL(3) DIMENSION CURR(512),IPK(4) REAL*4 NK INTEGER NTERMS INTEGER IDSRN DIMENSION IANS(2) DIMENSION LX(512),Y(512),WEIGHT(512),A(50),SIGMAA(50) INTEGER*4 ILOWER(4,8),IUPPER(4,8),IFIT(4) DIMENSION NCOL(4) COMMON/SCID/VO,JNE character*4 TITLE(16),TIT(2) character*4 VOY(2) character*4 PLNT1(2) character*4 PLNT2(2) character*4 PLNT3(2) character*4 PLNT4(2) character*8 CUP(4) character*8 DATE,TIME C REAL*4 voo(3) COMMON/SPECIE/AM(15),Z(15) INTEGER*2 MASSF(15),NQF(15) DIMENSION DEN(15),V(3,15),WPAR(15),WPER(15) DIMENSION DDEN(15),DV(3,15),DWPAR(15),DWPER(15) character*8 SFIT(10) LOGICAL*4 LFIRST REAL*4 A4(2) C REAL*4 ZFFFF/ZFFFF/ EQUIVALENCE (DATE,A4(1)) data ICPU/0/,JCPU/0/ DATA CUPN/-.2962,-.1710,-.9397,+.2962,-.1710,-.9397, * +0.0,.3420,-.9397,.7310,.6816,-.0349/ data BARR /645*0./ DATA TIT(2)/'BIMA'/,TIT(1)/' MA'/ DATA TITLE(2)/'XWEL'/,TITLE(3)/'LIAN'/ DATA TITLE(4)/' SIM'/,TITLE(5)/'ULAT'/,TITLE(6)/'ION '/ DATA TITLE(7)/'OF V'/ DATA TITLE(11)/' ON '/ data VOY /'1 AT','2 AT'/ data PLNT1 /' JUP','ITER'/,PLNT2 /' SAT','URN '/ data PLNT3 /' UR','ANUS'/,PLNT4 /' NEP','TUNE'/ data CUP /'A-CUP: ','B-CUP: ','C-CUP: ','D-CUP:'/ data SFIT /'PMODEL ','XJC ','XJC(UR) ','ABSO ', +'ABSOM ','BILAB ','COLD APX','XCUR ','CUPEXP ', +'XTEST '/ data LFIRST/.TRUE./ data voo /3*0.0/ c DATA PRESET/Z7FFFFFFF/ DATA PRESET/428244919/ C INTEGER MASK/ZFF40FFFF/ C ******************************************************************** C * * C * KDSRN( 6)=0 TURNS OFF DETAIL IN XMJSFT * C * KDSRN( 7)=0 TURNS OFF DETAIL IN XANAL AND "ANSWER" DETAIL * C * KDSRN(10)=0 TURNS OFF "I ( ) CALCULATED" MESSAGES * C * KDSRN(11)=0 TURNS OFF "PLSBEG" DETAIL * C * KDSRN(12)=0 TURNS OFF "VGRLOG" DETAIL * C * KDSRN(14)=0 TURNS OFF DIAGNOSTICS IN "XFNCDV" * C * * C ******************************************************************** C C RES ("RESULT") IS THE STANDARD DATA OUTPUT ARRAY NORMALLY CALLED C ANS ("ANSWER") C C C ******************************************************************** C * START EXECUTION OF CODE HERE * C ******************************************************************** C C C IANS(1) NUMBER OF WORDS IN FILLED "ANSWER" ARRAY C IANS(2) TALLY OF NUMBER OF CALLS TO "PLSANL" (0 IN THIS USAGE) C IPK(I) CHANNEL NUMBER OF PEAK CURRENT IN I-TH CUP C C IANS(1) AND IANS(2) ARE INTEGER*2 WORDS; THEY CONTAIN THE SAME C INFORMATION AS THE REAL*4 WORD RES(1) C C IPK(1) THROUGH IPK(4) ARE INTEGER*2 WORDS; THEY CONTAIN THE SAME C INFORMATION AS THE REAL*4 WORDS RES(2) AND RES(3) C CALL TIMER IPRT=6 CALL MVC(RES(1),IANS(1),4) CALL MVC(RES(2),IPK(1),8) C C IANS(1) SET TO 4 IN KNTCUR C IANS(1) SET TO 15 IN GETFLD C C CLEAR UPPER PART OF "RES" ARRAY C IF (.NOT. LFIRST) GO TO 90 DO 100 K=26,129 RES(K)=PRESET 100 CONTINUE DO 90 K=133,150 RES(K)=PRESET 90 CONTINUE C C INCREASE IANS(1) TO 150 C IANS(1)=150 CALL MVC(IANS(1),RES(1),4) C LQUAL = .FALSE. C C RMARK(1) IS THRESHOLD CURRENT C XCRIT=AMAX1(CSIG*CNOISE,RMARK(1)) C C B IS IN S/C COORDINATES HERE; ROTATED TO CYLINDRICAL IN XPARM C CALL MVC(RES(5),B(1),12) C C ********************************************************************* C * READ IN CONTROL INFORMATION FOR ANALYSIS * C ********************************************************************* C IOVER=0 94 WRITE (6,999) 999 FORMAT (1H1,5X,'** PARAMETERS FOR USE IN ANALYSIS **'//) IF (LFIRST) GO TO 1995 c 600 REWIND 5 600 continue WRITE (6,1999) 1999 FORMAT (1H ,'USE PREVIOUS FIT PARAMETERS AGAIN? 1=YES, 0=NO') READ (5,*,END=600,Err=600) IAG IF (IAG .EQ. 0) GO TO 610 IF (IAG .EQ. 1) GO TO 815 GO TO 600 C C READ LOGICAL FLAGS FOR FIT C c 610 REWIND 5 610 continue 1995 WRITE (6,800) 800 FORMAT (1H ,'INPUT NUMBER OF COMPONENTS TO FIT TO SPECTRA'//) READ (5,*,END=610,Err=610) NCOMP IF (NCOMP .LE. 15) GO TO 1802 WRITE (6,1800) 1800 FORMAT (1H ,//'** THE NUMBER OF COMPONENTS CANNOT EXCEED 15 **'//) GO TO 1995 1802 DO 775 J=1,NCOMP LN(J)=.TRUE. 775 CONTINUE c 620 REWIND 5 620 continue WRITE(6,770) 770 FORMAT(1H ,'DO YOU WISH TO HAVE DENSITY AS A FIT PARAMETER?', * ' (T OR F)') READ(5,*,END=620,Err=620) LDEN IF(LDEN) GO TO 789 DO 780 J=1,NCOMP LN(J)=.FALSE. 780 CONTINUE GO TO 790 789 IF(NCOMP.NE.1)GO TO 788 LN(NCOMP)=.TRUE. GO TO 790 788 WRITE(6,747)NCOMP 747 FORMAT(1H ,'INDICATE WHICH DENSITIES TO FIT (T OR F FOR EACH OF ', * 'THE ',I1,' SPECIES)') READ(5,*,END=7000)(LN(N),N=1,NCOMP) 790 CONTINUE DO 802 M=1,NCOMP WRITE (6,803) 803 FORMAT (1H ,'INPUT LOGICAL VARIABLES (T OR F ) TO SELECT FIT PARAM +ETERS'/) c 630 REWIND 5 630 continue WRITE (6,804) M 804 FORMAT (1H ,'FIT VR, VPHI, VZ FOR COMPONENT ',I2,' ?') READ (5,*,END=630) LV(1,M),LV(2,M),LV(3,M) c 640 REWIND 5 640 continue WRITE (6,805) 805 FORMAT (1H ,'FIT THERMAL SPEED FOR THIS COMPONENT ?') READ (5,*,END=640) LW(M) C C READ MASSES AND CHARGES OF COMPONENTS C c 650 REWIND 5 650 continue WRITE (6,807) M 807 FORMAT (1H ,'INPUT MASS AND CHARGE STATE ( INTEGERS ) OF COMPONENT + ',I2) READ (5,*,END=650) MASSF(M),NQF(M) AM(M)=MASSF(M) Z(M)=NQF(M) 802 CONTINUE C C READ CONTROL PARAMETERS FOR FIT C WRITE (6,810) 810 FORMAT (1H //) c 660 REWIND 5 660 continue WRITE (6,808) 808 FORMAT (1H ,'USE ISOTROPIC MAXWELLIANS IN FIT ?') READ (5,*,END=660) LM LTH=.TRUE. IF (NCOMP .EQ. 1) GO TO 1994 DO 2671 M=1,NCOMP IF ( .NOT. LW(M)) GO TO 1671 2671 CONTINUE GO TO 1994 c 1671 REWIND 5 1671 continue WRITE (6,809) 809 FORMAT (1H ,'USE COMMON TEMPERATURE (ALL THERMAL SPEEDS ARE NOT', +' FIT PARAMETERS) ?') READ (5,*,END=1671) LTH 1994 FLAM=.100 ALAM=.001 IF (IOVER .EQ. 0) GO TO 815 WRITE (6,811) 811 FORMAT (1H ,'MARQUARDTS PARAMETER (FLAM) IS SET TO 0.100') WRITE (6,812) 812 FORMAT (1H ,'THE PARTIAL DERIVATIVE PARAMETER (ALAM) IS SET TO 0.0 +01') c 680 REWIND 5 680 continue WRITE (6,813) 813 FORMAT (1H ,'OVERRIDE THESE VALUES ? 1=YES, 0=NO ') READ (5,*,END=680) IOVER IF (IOVER .EQ. 0) GO TO 815 IF (IOVER .EQ. 1) GO TO 690 GO TO 680 c 690 REWIND 5 690 continue WRITE (6,816) 816 FORMAT (1H ,'INPUT FLAM ') READ (5,*,END=690) FLAM c 700 REWIND 5 700 continue WRITE (6,817) 817 FORMAT (1H ,'INPUT ALAM ') READ (5,*,END=700) ALAM 815 CONTINUE C C WRITE OUT FIT PARAMETERS AND OPTIONS C C IF (IAG .EQ. 1) NCOMP=RES(48) WRITE (6,900) 900 FORMAT(1H ,'NCOMP',1X,'CHARGE',2X,'MASS',2X,'FIT N',1X, +'FIT VR',1X,'FIT VPHI',1X,'FIT VZ',1X,'FIT W'/) DO 910 M=1,NCOMP 910 WRITE (6,920) M,NQF(M),MASSF(M),LN(M),(LV(M1,M),M1=1,3),LW(M) 920 FORMAT (1X,I5,1X,2I6,L6,L7,L9,L7,L6,5X,2I5) WRITE (6,901) 901 FORMAT (1H /) WRITE (6,921) FLAM,ALAM 921 FORMAT (1H ,'FLAM = ',F10.3,' ALAM = ',F10.3) WRITE (6,901) WRITE (6,922) LM 922 FORMAT (1H ,'ISOTROPIC MAXWELLIAN ',L5) IF (NCOMP .NE. 1) WRITE (6,923) LTH 923 FORMAT (1H ,'SAME TEMPERATURE ',L5) WRITE (6,901) c 710 REWIND 5 710 continue WRITE (6,924) 924 FORMAT (1H ,'CHANGE FIT OPTIONS? 1=YES, 0=NO') READ (5,*,END=710) IOVER IF (IOVER .EQ. 0) GO TO 720 IF (IOVER .EQ. 1) GO TO 94 GO TO 710 C C FILL LX AND Y WITH CHANNEL NUMBERS AND Y WITH CORRESPONDING C CURRENTS WHICH ARE TO BE USED IN THE FIT C C C ****************************************************** C * MAIN LOOP OVER ALL 4 CUPS * C ****************************************************** C 720 CALL CLRTEK IF (LFIRST) GO TO 1996 c 730 REWIND 5 730 continue WRITE(6,1998) 1998 FORMAT (1H ,'FIT SAME CHANNELS AS BEFORE? 1=YES, 0=N0') READ(5,*,END=730) IAG IF (IAG .EQ. 0) GO TO 1996 IF (IAG .EQ. 1) GO TO 106 GO TO 730 1996 IADD=0 JADD=0 DO 102 J=1,4 c 740 REWIND 5 740 continue WRITE (6,5990) J 5990 FORMAT (1H ,'FIT DATA IN CUP',I2,' ? 1=YES, 0=NO') READ (5,*,END=740) IFIT(J) IF (IFIT(J) .EQ. 0) GO TO 104 IF (IFIT(J) .EQ. 1) GO TO 5991 GO TO 740 5991 CONTINUE LOWER=0 JX=(J-1)*NE IDSRN=KDSRN(7) IF (IDSRN .NE.0) WRITE (7,6000) 6000 FORMAT (1H0,'FROM XANAL',2X,'IADD',5X,'L',10X,'Y',5X,'WEIGHT',5X, +'WT2'/) c 750 REWIND 5 750 continue WRITE(6,430) 430 FORMAT(1H ,'HOW MANY CONTIGUOUS SETS OF CHANNELS?') READ(5,*,END=750) JSETS(J) JSET=JSETS(J) DO 103 I=1,JSET c 760 REWIND 5 760 continue WRITE (6,5992) I 5992 FORMAT (1H ,'INPUT LOWER AND UPPER CHANNELS FOR SET ',I1) READ (5,*,END=760) ILOWER(J,I),IUPPER(J,I) L1=MAX0(ILOWER(J,I)*1,LOWER*1+1) L2=MIN0(IUPPER(J,I)*1,NE*1) IF(L1.GT.NE) L1=NE IF(L2.LT.L1) L2=L1 DO 200 L=L1,L2 RCURR=CURR(JX+L) IF(RCURR.LE.RMARK(1)) GO TO 200 C C RMARK(1) IS ZERO DN LEVEL C IF(RCURR.EQ.RMARK(3)) GO TO 200 C C RMARK(3) IS SATURATION CURRENT C IADD=IADD+1 LX(IADD)=L C C NK(J,MD) CONVERTS TO UNITS OF KM/(SEC CM**3) C ********************** N. B. ********************************* C * THE NK MATRIX HAS BEEN THE SOURCE OF INFINITE PROBLEMS; * C * SINCE SIMULATED CURRENTS ARE BEING GENERATED WITH PMODEL * C * IN FEMTOAMPS AND RCURR IS IN FEMTOAMPS, WE WILL ALSO PUT * C * THE Y MATRIX IN FEMTOAMPS. * C ************************************************************** C Y(IADD)=RCURR ICHAN=0 C IF (LWT) ICHAN=L C C ********************************************************************* C * * C * *** STATISTICAL WEIGHTS FOR ERROR ANALYSIS *** * C * * C * THE THERMAL NOISE (IN FEMTOAMPS) IS CNOISE WHERE * C * CNOISE=FAC(ICLK)*SIGNOI * C * SIGNOI=35.0 FEMTOAMPS * C * FAC(ICLK)=SQRT(CLOCK(3)/CLOCK(ICLK)) * C * CLOCK(ICLK) IS THE INTEGRATION TIME IN SECONDS * C * * C * ICLK CLOCK(ICLK) FAC(ICLK) * C * 1 0.03 5.6 * C * 2 0.21 2.1 * C * 3 0.93 1.0 * C * * C * FOR GS-3 ICLK=2 AND CNOISE=73.5 FEMTOAMPS * C * THE DIGITIZATION ERROR IS SQRT(FNOIS)*RCURR IN FEMTOAMPS WHERE * C * FNOIS=(10**1/64-1)**2/12 = 1.118E-4 IS THE APPROPRIATE * C * FACTOR * C * * C ********************************************************************* C C WEIGHT(IADD) ONLY INCLUDES THE THERMAL NOISE AND IS THE SAME C FOR ALL MEASURED CURRENTS C WEIGHT(IADD)=1./(CNOISE*CNOISE) C C WT2(IADD) INCLUDES BOTH THERMAL NOISE AND DIGITIZATION ERROR C WT2(IADD)=1./(CNOISE*CNOISE+FNOIS*RCURR*RCURR) IF (IDSRN .EQ. 0) GO TO 200 WRITE (7,6100) IADD,L,Y(IADD),WEIGHT(IADD),WT2(IADD) 6100 FORMAT (1H ,11X,2I6,3E11.3) 200 CONTINUE 103 CONTINUE 104 CONTINUE NCOL(J)=IADD-JADD JADD=IADD 102 CONTINUE GO TO 162 106 CONTINUE IADD=0 JADD=0 LOWER=0 DO 357 J=1,4 IF(IFIT(J).EQ.0)GO TO 366 JX=(J-1)*NE JSET=JSETS(J) DO 367 I=1,JSET L1=MAX0(ILOWER(J,I),LOWER+1) L2=MIN0(IUPPER(J,I)*1,NE*1) IF(L1.GT.NE)L1=NE IF(L2.LT.L1)L2=L1 DO 377 L=L1,L2 RCURR=CURR(JX+L) IF(RCURR.LE.RMARK(1))GO TO 377 IF(RCURR.EQ.RMARK(3))GO TO 377 IADD=IADD+1 LX(IADD)=L Y(IADD)=RCURR WEIGHT(IADD)=1./(CNOISE*CNOISE) WT2(IADD)=1./(CNOISE*CNOISE+FNOIS*RCURR*RCURR) 377 CONTINUE 367 CONTINUE 366 CONTINUE NCOL(J)=IADD-JADD JADD=IADD 357 CONTINUE 162 CONTINUE C C ******************************************************* C * END MAIN LOOP OVER CUPS * C ******************************************************* C C C INITIALIZE FIT PARAMETERS C CALL XPARM(NCOMP,NQF,MASSF,CURR,DEN,V,WPAR,WPER,RES) C C PUT INITIAL VALUES INTO "A" ARRAY C 5555 CONTINUE CALL XPACK(NCOMP,DEN,V,WPAR,WPER,A,NTERMS) CALL CLRTEK IF (KDSRN(6) .NE. 0) WRITE (6,2005) NTERMS 2005 FORMAT(1H ,'NTERMS = ',I10) FLAMD=FLAM ICHI=0 C C CALL FITTING ROUTINE C MODE=MD IDSRN=KDSRN(6) IF(MD.EQ.1) CALL XMJSFT(LX,Y,WEIGHT,NCOL,NTERMS,NCOMP,ICHI,A, * SIGMAA,FLAMD,CHISQR,VJWL,ICALL,IQUAL,IDSRN,EPS,IEND,NE) IF (MD .EQ. 1 .AND. ICPLT .EQ. 5) GO TO 404 IF(MD.EQ.2) CALL XMJSFT(LX,Y,WEIGHT,NCOL,NTERMS,NCOMP,ICHI,A, * SIGMAA,FLAMD,CHISQR,VJWM,ICALL,IQUAL,IDSRN,EPS,IEND,NE) GO TO 401 C C CAN DO M-MODES ONLY WITH ABSOM C 404 WRITE (6,405) 405 FORMAT (1H ,'L-MODE SPECTRA CANNOT BE ANALYZED WITH ABSOM') 401 IDSRN=KDSRN(7) IF ( CHISQR.GT.0. ) RES(82) = ALOG10( CHISQR) IF (CH2 .GT. 0.) RES(83)=ALOG10(CH2) RES(84) = ICALL RES(85) = IQUAL DO 107 JCUP=1,4 107 RES(85+JCUP)=NCOL(JCUP) C C IQUAL=0 - FIT CRITERIA SATISFIED C IQUAL=1 - FIT CRITERIA NOT SATISFIED C IQUAL=2 - FIT FAILED IN XFNCDV OR XMJSFT C IQUAL=3 - NO POINTS GIVEN TO XMJSFT TO FIT C IQUAL=4 - FIT CRITERIA NOT SATISFIED; CHISQR INCREASING ON LAST C ITERATION C IF (IQUAL .EQ. 1) WRITE(6,421) IF (IQUAL .EQ. 2) WRITE(6,422) IF (IQUAL .EQ. 3) GO TO 3005 421 FORMAT(1H ,'IQUAL=1 - FIT CRITERIA NOT SATISFIED') 422 FORMAT(1H ,'IQUAL=2 - FIT FAILED IN XFNCDV OR XMJSFT') c 771 REWIND 5 771 continue WRITE(6,345) 345 FORMAT(1H ,'OUTPUT PARTICLE VELOCITIES IN', * ' S/C COORDINATES? 1=YES, 0=NO') READ(5,*,END=771) IDOO IF (IDOO .NE. 0 .AND. IDOO .NE. 1) GO TO 771 c 772 REWIND 5 772 continue WRITE(6,643) 643 FORMAT(1H ,'OUTPUT PARTICLE VELOCITIES IN', * ' RTN COORDINATES? 1=YES, 0=NO') READ(5,*,END=772) JDOO CALL UBELL CALL UPAUSE CALL CLRTEK IF (JDOO .NE. 0 .AND. JDOO .NE. 1) GO TO 772 C C UNPACK FIT PARAMETER ARRAY AND UNCERTAINTIES C CALL CLRTEK CALL XUNPCK(A,SIGMAA,NCOMP,DEN,DDEN,V,DV,WPAR,DWPAR +,WPER,DWPER) C C **************************************************************** C * * C * DDEN ETC. CONTAIN ONLY THE APPROPRIATE VALUES OF THE * C * SIGMAA ARRAY. TO ESTIMATE THE TRUE UNCERTAINTIES IN * C * THE FIT PARAMETERS, THESE VALUES SHOULD BE MULTIPLIED * C * BY THE SQUARE ROOT OF THE REDUCED CHI SQUARE OF THE * C * FIT. THIS TECHNIQUE ESTIMATES THE ERROR IN THE DATA * C * POINTS BY USING THE VARIANCE IN THE FIT (SEE BEVINGTON). * C * * C **************************************************************** C IF(IQUAL.EQ.0) WRITE(6,676) 676 FORMAT(1H0,53X,'*** THIS IS A BEST FIT ***') CALL DATI(DATE,TIME) WRITE(6,32)TIME,DATE 32 FORMAT(1H0,'FIT COMPLETED AT ',A8,' ON ',A8) CALL TIMER CALL HOTIME C C DETERMINE COST OF SIMULATION C JCPU=ITT-ICPU ICPU=ITT FTM=JCPU*26.04166E-6/60. WRITE (6,2798)FTM 2798 FORMAT(1H ,'CPU TIME USED = ',F10.5,2X,'MINUTES'/) DFTM=FTM*24. WRITE (6,2797) DFTM 2797 FORMAT (1H ,'THIS FIT COST x ',F7.2,' (PRIME TIME RATE)') C C WRITE OUT SPACECRAFT AND PLANET IDENTIFICATION C IF(.NOT.LM) TITLE(1)=TIT(2) TITLE(1)=TIT(1) C IF(VO .EQ. 453.) TITLE(8)=VOY(1) IF(VO .EQ. 372.) TITLE(8)=VOY(2) GO TO(25,30,35,40),NUMFOR 40 TITLE(9)=PLNT4(1) TITLE(10)=PLNT4(2) GO TO 45 35 TITLE(9)=PLNT3(1) TITLE(10)=PLNT3(2) GO TO 45 30 TITLE(9)=PLNT2(1) TITLE(10)=PLNT2(2) GO TO 45 25 TITLE(9)=PLNT1(1) TITLE(10)=PLNT1(2) 45 CONTINUE CALL MVC(TIMESR(1),TITLE(12),20) WRITE(6,82)(TITLE(I),I=1,16) 82 FORMAT(1H0,'INFORMATION FROM ',16A4) WRITE(6,83) 83 FORMAT(//) WRITE(6,84)SFIT(ICPLT) 84 FORMAT(1X,'CUP CHANNELS USED IN FIT BY ',A8) C C DO LOOP OVER CUPS C DO 50 I=1,4 IF(IFIT(I).NE.0) GO TO 51 WRITE(6,86)CUP(I) GO TO 50 51 CONTINUE C C DO LOOP OVER SPECIES C JSET=JSETS(I) DO 60 J=1,JSET IF(J.NE.1) GO TO 70 WRITE(6,80)CUP(I),ILOWER(I,J),IUPPER(I,J) GO TO 60 70 WRITE(6,81)ILOWER(I,J),IUPPER(I,J) 80 FORMAT(1H0,A8,'CHANNELS ',I3,' -',I3) 81 FORMAT(1H ,8X,'CHANNELS ',I3,' -',I3) 60 CONTINUE 86 FORMAT(1H0,A8,'*** NO FIT PERFORMED ***') 50 CONTINUE WRITE(6,85)CHISQR 85 FORMAT(1H0,'CHI SQUARE = ',E11.3) C SQCHI=SQRT(CHISQR) WRITE (6,601) 601 FORMAT (1H /) WRITE (6,602) TIMESR 602 FORMAT (1H ,'PARAMETERS AND 1-SIGMA UNCERTAINTIES FROM FIT TO SP', 1'ECTRUM AT ',5A4/) WRITE(6,599) SFIT(ICPLT) 599 FORMAT(1H ,'CURRENTS WERE GENERATED USING ',A8) IF (LD) WRITE (6,2793) 2793 FORMAT (1H ,'DOUBLE PRECISION CODE WAS USED'/) IF ( .NOT. LD) WRITE (6,2794) 2794 FORMAT (1H ,'SINGLE PRECISION CODE WAS USED'/) IF (LM) WRITE (6,603) 603 FORMAT (1H ,' MASS NQ',7X,'DEN',6X,'XDEN',8X,'VR',7X,'XVR', 16X,'VPHI',5X,'XVPHI',8X,'VZ',7X,'XVZ',9X,'W',8X,'XW'/) IF ( .NOT. LM) WRITE (6,1603) 1603 FORMAT (1H ,' MASS NQ',7X,'DEN',6X,'XDEN',8X,'VR',7X,'XVR', 16X,'VPHI',5X,'XVPHI',8X,'VZ',7X,'XVZ',6X,'WPAR',5X,'XWPAR', 26X,'WPER',5X,'XWPER'/) DO 605 M=1,NCOMP XDEN(M)=SQCHI*DDEN(M) XVR(M)=SQCHI*DV(1,M) XVPHI(M)=SQCHI*DV(2,M) XVZ(M)=SQCHI*DV(3,M) XWPAR(M)=SQCHI*DWPAR(M) XWPER(M)=SQCHI*DWPER(M) IF (LM) WRITE (6,604) MASSF(M),NQF(M),DEN(M),XDEN(M),V(1,M), 1XVR(M),V(2,M),XVPHI(M),V(3,M),XVZ(M),WPAR(M),XWPAR(M) 604 FORMAT (1H ,2I5,10F10.4) IF ( .NOT. LM) WRITE (6,1604) MASSF(M),NQF(M),DEN(M),XDEN(M), 1V(1,M),XVR(M),V(2,M),XVPHI(M),V(3,M),XVZ(M), 2WPAR(M),XWPAR(M),WPER(M),XWPER(M) 1604 FORMAT (1H ,2I5,12F10.4) 605 CONTINUE C C----- SBLK(40-42) IS S/C VEL. IN S/C COOR. C DO 439 J=1,NCOMP CALL VTRNS(SWVSV,TCYSSC,V(1,J),voo) CALL VTRNS(VSWRTN,TSCRTN,SWVSV(1),voo) DO 333 I=1,3 SWV(I,J)=SWVSV(I)-SBLK(39+I) VRTN(I,J)=VSWRTN(I) 333 CONTINUE 439 CONTINUE IF(IDOO.EQ.0)GO TO 765 WRITE(6,321) 321 FORMAT(1H0,'PARTICLE VELOCITIES IN S/C COORDINATES (KM/SEC)') WRITE(6,322) 322 FORMAT(1H0,5X,'NCOMP',5X,'VX',8X,'VY',8X,'VZ',/) DO 531 J=1,NCOMP WRITE(6,334)J,SWV(1,J),SWV(2,J),SWV(3,J) 334 FORMAT(8X,I1,5X,F7.2,3X,F7.2,3X,F7.2) 531 CONTINUE 765 CONTINUE IF(JDOO.EQ.0)GO TO 766 WRITE(6,323) 323 FORMAT(1H0,'PARTICLE VELOCITIES IN RTN COORDINATES (KM/SEC)', * ' *** CORRECTED FOR S/C MOTION ***') WRITE(6,324) 324 FORMAT(1H0,5X,'NCOMP',5X,'VR',8X,'VT',8X,'VN',/) DO 753 J=1,NCOMP WRITE(6,343)J,VRTN(1,J),VRTN(2,J),VRTN(3,J) 343 FORMAT(8X,I1,5X,F7.2,3X,F7.2,3X,F7.2) 753 CONTINUE 766 CONTINUE 346 CONTINUE WRITE(6,666) 666 FORMAT(1H0,'PROJECTION OF FITTED VELOCITIES ALONG CUP NORMALS', * ' (KM/SEC)') WRITE(6,667) 667 FORMAT(1H0,5X,'NCOMP',5X,'VA',8X,'VB',8X,'VC',8X,'VD',/) DO 670 J=1,NCOMP DO 671 I=1,4 VCUPN(I,J)=-(SWV(1,J)*CUPN(1,I)+SWV(2,J)*CUPN(2,I)+ * SWV(3,J)*CUPN(3,I)) 671 CONTINUE WRITE(6,672)J,(VCUPN(I,J),I=1,4) 672 FORMAT(8X,I2,5X,F7.2,3X,F7.2,3X,F7.2,3X,F7.2) 670 CONTINUE C C OPTION TO PLOT SIMULATED CURRENTS (TOTAL FROM ALL SPECIES) C c 781 REWIND 5 781 continue WRITE(6,370) 370 FORMAT(1H0,'PLOT FIT CURRENTS FROM BEST FIT ? 1=YES, 0=NO') READ(5,*,END=781) IPLT IF (IPLT .EQ. 0) GO TO 360 IF (IPLT .GE. 1 .OR. IPLT .LE. 4) GO TO 791 GO TO 781 791 CONTINUE IF (IQUAL .EQ. 0 .OR. IQUAL .EQ. 1) GO TO 792 IF (IQUAL .EQ. 4) GO TO 793 WRITE (6,794) 794 FORMAT (1H ,'FAILURE MODE PRECLUDES PLOTTING. TRY NEW INPUT ', *'PARAMETERS.') GO TO 360 793 CONTINUE WRITE (6,796) 796 FORMAT (1H ,'CHISQR INCREASED ON THE LAST ITERATION.'/' TO PLOT ', *'THE BEST FIT CURRENTS AT LEAST ONE MORE ITERATION IS NECESSARY') c 1797 REWIND 5 1797 continue WRITE (6,797) 797 FORMAT (1H ,'DO YOU STILL WISH TO PLOT THE BEST FIT? 1=YES, 0=NO') READ (5,*,END=1797) ICON IF (ICON .EQ. 0) GO TO 360 IF (ICON .EQ. 1) GO TO 5555 GO TO 1797 792 NCMP=NCOMP WRITE(6,877)NCMP 877 FORMAT(1H ,2X,'NCMP = ',I2) ME=0 MP=MD-1 C C ZERO BARR C BARR(1,1)=0.0 CALL MVC(BARR(1,1),BARR(2,1),2576) DO 500 JCUP=1,4 DO 510 L=1,NE DO 520 N=1,NCMP BARR(L,JCUP)=BARR(L,JCUP)+BTOT(L,JCUP,1,N) 520 CONTINUE 510 CONTINUE 500 CONTINUE CALL veccpy (VEL(1),V(1,1)) CALL PLTSIM(BARR,MP,ME,NCMP,AM,Z,VEL,B,DEN,WPAR,WPER) WRITE(6,749) 749 FORMAT(' ') 360 CONTINUE C C THIS SECTION OUTPUTS THE FIT PARAMETERS TO THE USERS A-DISK C ON FILE FT40F001 (IF FILE IS NOT EMPTY THE NEW FIT C PARAMETERS ARE ADDED TO THE END OF THE FILE) C C WRITE(6,851) C 851 FORMAT(1H ,'OUTPUT THIS FIT TO YOUR DISK? 1-9 = YES, 0,CR = NO') C READ(5,*,END=860) KEEP C IF(KEEP.EQ.0) GO TO 860 C WRITE(6,9091) C9091 FORMAT(1H ,'INPUT YOUR INITIALS') C READ(5,9092)UINIT C9092 FORMAT(A4) C PUT SPECTRUM TIME IN PROPER FORM C TIMESR(2)=SHFTL(TIMESR(2),8) C TIMESR(4)=AND(TIMESR(4),MASK) C DECIDES IF L OR M-MODE C IF(MD.EQ.2)MMLL=LLMM(1) C IF(MD.EQ.1)MMLL=LLMM(2) C B FIELD AND CHISQR ARE SPLIT INTO A NUMBER AND C EXPONENT TO SAVE SPACE IN OUTPUT FILE C IE1=0 C IE2=0 C IE3=0 C IF(ABS(B(1)).GE.10)IE1=1 C IF(ABS(B(2)).GE.10)IE2=1 C IF(ABS(B(3)).GE.10)IE3=1 C IF(ABS(B(1)).GE.100)IE1=2 C IF(ABS(B(2)).GE.100)IE2=2 C IF(ABS(B(3)).GE.100)IE3=2 C IF(ABS(B(1)).GE.1000)IE1=3 C IF(ABS(B(2)).GE.1000)IE2=3 C IF(ABS(B(3)).GE.1000)IE3=3 C IC1=0 C IF(CHISQR.GT.1.E1)IC1=1 C IF(CHISQR.GT.1.E2)IC1=2 C IF(CHISQR.GT.1.E3)IC1=3 C IF(CHISQR.GT.1.E4)IC1=4 C IF(CHISQR.GT.1.E5)IC1=5 C IF(CHISQR.GT.1.E6)IC1=6 C IF(CHISQR.GT.1.E7)IC1=7 C IF(CHISQR.GT.1.E8)IC1=8 C IF(CHISQR.GT.1.E9)IC1=9 C CHISQU=CHISQR/10**IC1 C B1=B(1)/10**IE1 C B2=B(2)/10**IE2 C B3=B(3)/10**IE3 C PUTS FIT DATE IN PROPER FORM (REMOVES SLASHES) C AM1=SHFTL(A4(1),24) C AD=AND(A4(2),ZFFFF) C AD=SHFTL(A4(2),16) C WRITE FIRST LINE IN FILE(TRAJ INFO,TIMES,B,ETC.) C WRITE(40,852) SSSS,TIMESR,SBLK(1),SBLK(19),B1,EEEE,IE1,B2,EEEE,IE2 C 1,B3,EEEE,IE3,MMLL,A4(1),AM1,A4(2),AD,UINIT, C 2ICPLT,CHISQU,EEEE,IC1,NCOMP C 852 FORMAT(A1,A4,A3,3A4,2F7.3,F5.2,A1,I1,F5.2,A1,I1,F5.2,A1,I1, C 1A2,1X,A2,A1,A1,A2,1X,A3,I2,F5.3,A1,I1,I2) C WRITE CHANNELS FIT IN EACH CUP TO FILE C JASET=JSETS(1) C JBSET=JSETS(2) C JCSET=JSETS(3) C JDSET=JSETS(4) C IF(IFIT(1).EQ.1)WRITE(40,853)CUP(1),JSETS(1),((ILOWER(1,I) C 1,IUPPER(1,I)),I=1,JASET) C IF(IFIT(1).EQ.0)WRITE(40,853)CUP(1),IFIT(1) C IF(IFIT(2).EQ.1)WRITE(40,853)CUP(2),JSETS(2),((ILOWER(2,I) C 1,IUPPER(2,I)),I=1,JBSET) C IF(IFIT(2).EQ.0)WRITE(40,853)CUP(2),IFIT(2) C IF(IFIT(3).EQ.1)WRITE(40,853)CUP(3),JSETS(3),((ILOWER(3,I) C 1,IUPPER(3,I)),I=1,JCSET) C IF(IFIT(3).EQ.0)WRITE(40,853)CUP(3),IFIT(3) C IF(IFIT(4).EQ.1)WRITE(40,853)CUP(4),JSETS(4),((ILOWER(4,I) C 1,IUPPER(4,I)),I=1,JDSET) C IF(IFIT(4).EQ.0)WRITE(40,853)CUP(4),IFIT(4) C WRITE FIT PARAMETERS FOR EACH ION SPECIES INTO FILE C DO 855 JI=1,NCOMP C WRITE(40,856)FFFF,JI,MASSF(JI),NQF(JI),DEN(JI),XDEN(JI),V(1,JI), C 1XVR(JI),V(2,JI),XVPHI(JI) C 856 FORMAT(A1,I1,I3,I2,6E10.3) C 855 WRITE(40,857)GGGG,JI,V(3,JI),XVZ(JI),WPAR(JI),XWPAR(JI),WPER(JI), C 1XWPER(JI) C 857 FORMAT(A1,I1,6E10.3) C 853 FORMAT(A1,21I3) C C INSERT PHYSICAL PARAMETERS AND THEIR UNCERTAINTIES INTO RES C (THE "ANSWER ARRAY") C 860 RES(48)=NCOMP IWRS=48 DO 4000 N=1,NCOMP II=(N-1)*12+IWRS RES(II+1)=DEN(N) RES(II+2)=DDEN(N) CALL MVC(V(1,N),RES(II+3),12) CALL MVC(DV(1,N),RES(II+6),12) RES(II+9)=WPAR(N) RES(II+10)=WPER(N) RES(II+11)=DWPAR(N) RES(II+12)=DWPER(N) 4000 CONTINUE 7000 CONTINUE LFIRST=.FALSE. C C OUTPUT C IDSRN= KDSRN(7) IF (IDSRN .EQ. 0) RETURN WRITE(6,3003) (LPLS(J),J=9,16),LANIS,LALPHI,NQ,MASS,KEY,NCOM,NVAR 3003 FORMAT(1H0,3(4L1,1X),L2,2(4I1,2X),5X,4I3,5X,2I3,5X,2I3) WRITE(6,3004) IANS,IPK,XCRIT 3004 FORMAT(1H ,2I10,'/',4I5,'/',E12.3) DO 3002 J=1,15 J1=(J-1)*10+1 J2=J1+9 WRITE(6,3001) (N,RES(N),N=J1,J2) 3001 FORMAT(1H ,10(I3,F9.2)) 3002 CONTINUE RETURN 3005 CONTINUE IF (IDSRN .NE. 0) WRITE (6,3006) IQUAL,XCRIT 3006 FORMAT (1H0,'FIT NOT DONE. IQUAL=',I4,' ; XCRIT=',E10.3) RETURN END