C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C ASFIT - ASYMMETRIC ROTOR FITTING PROGRAM, using Watson's reduced C Hamiltonian with terms up to decadic. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Principal features: C C - J UP TO JMAX=250 C - ALL CONSTANTS IN WATSON'S REDUCED HAMILTONIAN C UP TO AND INCLUDING DECADIC CAN BE FITTED FOR REDUCTIONS C 'A' AND 'S' AND REPRESENTATIONS 'I.r' AND 'III.r'. C - UP TO NLINES=5000 TRANSITIONS C - WEIGHTED FIT C - OPTIONAL CORRELATION COEFFICIENTS, MOST INLUENTIAL TRANSITIONS C AND CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO TRANSITION C FREQUENCIES (EVALUATED FROM DERIVATIVES) C - FIT WITH ZERO DEGREES OF FREEDOM (ie. N constants to N C frequencies) IS ALLOWED, BUT THE USER IS WARNED BY THE BELL C - OPTION TO FIT PEAK FREQUENCIES OF BLENDED LINES C - GENERATION OF QUANTUM NUMBERS OF COMPONENT TRANSITIONS FOR C MANY BAND TYPES C - ON-LINE DATA MODIFICATION (constants, transitions) DURING C THE FIT C - POSSIBILITY OF ANNOTATING EACH LINE WITH AN 8-BYTE COMMENT as C well as allowance for spacing/comment lines between C declared transitions C C This program, though much developed, owes much to the excellent ASFIT C of P.J.Mjoberg et al (ca 1972), which in turn seems to have carried C code written in Wilson's group. C C version 14.IV.2008 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- C C Modification history: C C 1979: port from IBM-360 to PDP-11 C 1980's: MERA-60 -> HP-150 -> PC's -> PRIME, SUN, SG C 1.06.95: allowance for annotations at ends of lines C 11.11.95: modification of data file creation from keyboard; C generation of oblate-type bands according to n=l/m generalization C 15.11.95: addition of deposition format for frequency table C 11.03.96: nomenclature I(+)(-),II(+)(-) and n=l/m scheme for II(-) C 12.03.96: modified line change option C 28.03.96: further flexibility in printout of derivatives C 10.04.96: explicit declaration of the handedness of representations used C 15.04.96: flagging of non-excluded lines with o-c > 3*exptal error C 11.04.97: input file='old', jmax=250, elimination of excl/annotations clash C 16.04.97: change in .ASF format, in output + debugging C 30.03.00: debugging of output C 7.04.01: overhaul of input/output C 10.01.02: addition of hjj/DFBETAS tests + major overhaul/debugging C 4.04.02: fitting of blended transitions C 4.01.03: various modifications to input and output C 7.04.03: .LIN output for SPFIT (preserving some comments for use by PIFORM) C 3.09.03: elimination of bug in fitting blended lines, line generation by C cloning + many incremental mods, including systematisation of C end-of-line comments and switch to 8-byte length C 12.08.04: modified line counting in data file, .PAR output + other mods C 3.12.04: modified output formats and data checks C 19.01.06: detection of potential split blends in dataset C 14.04.08: mods to SPFIT output C_____________________________________________________________________________ C C This program is standard F77 and compiles accordingly (tested with C quite a few FORTRANs, although not necessarily this particular version). C Any metacommands at the top of the listing beginning with a $ are for C MSFortran only and should be discarded on transfer to a different Fortran. C C Operating systems/compilers differ in handling of the old FORTRAN carriage C control, so some FORMAT statements may need to be modified - in particular C FORMATs number 801 and 802. C C Command for Compaq Visual Fortran compilation for all Pentium machines: C C df -arch=p5 -tune=p5 -fast -static -ccdefault:fortran asfit.for C C As above but smaler EXE: C C df -nopdbfile -nodebug -arch=p5 -tune=p5 -fast -static C -ccdefault:fortran asfit.for C C SG compilation: C C f77 -O2 -static -o asfit asfit.for C C_____________________________________________________________________________ C C DATA INPUT: C C A new data file is best generated by editing a previous one, which C should be self explanatory. The separating lines written by ASFIT into C the data set define the columns where the various numbers are to C be written (in right-justified form). C C The program provides its own limited editing facilities for C generating or modifying a data set. Generation from scratch is now rarely C used but the editing facility is still invaluable for making C modifications to the data as a result of a fit, and the program C allows full cyclic operation of C fitting->inspection of results->modification->fitting C C The following input conventions apply: C C - 1/0 convention for YES/NO (usually the answer is read in in fixed C format so an suffices for NO) C - constants are to be supplied in the form: 1,constant value C or: 0,constant value C depending on whether the constant is to be fitted or not C C - transitions are to be specified with the three rotational quantum C numbers J,K-1,K+1 of the upper state first, followed by the q.n's C of the lower state, the frequency, its error and fitting parameter: C 3,1,3,2,1,2,8669.6161,0.001,1 C the fitting parameter can be one of: C 0=do not fit C 1=fit C 2=fit blends, in which case intensity weighted mean frequency C of several calculated transitions is fitted to a single C measured frequency of the blended line. In this case C relative intensities of blended components follow C the fitting parameter, for example for unresolved C K2 doublet suspected of splitting greater than the frequency C error: C 3,2,2,2,2,1,8669.6161,0.05,2,1.0 C 3,2,1,2,2,0,8669.6161,0.05,2,1.0 C In the file output each component will have the same o-c C frequency, but also its individual relative intensity and C calculated frequency C C - termination of frequency input loop (including insert and C add modes) is by inputting 0's for every parameter eg. C 0,0,0,0,0,0,0,0 or ,,,,,,,, C C - exit from most modification options is by just pressing C the key C C If ASFIT is used to generate a data set from scratch (which may C sometimes be faster) use option 7 to save intermediate stages as well C as just before the first fitting cycle. C C * The number of transitions specified in the data file only serves as C the count of the number of lines at the moment of writing this file C by ASFIT. C On input the program reads up to the end of the data file, C IRRESPECTIVE of the number of transitions in the third line of file. C C * Empty lines can be placed between transition lines, and will C be translated as spaces in the same places in the printout C * Lines between transition lines beginning with C the '!' or '%' character are treated as comment lines - but only C those beginning with '!' will be echoed in the output C C_____________________________________________________________________________ C C The available 'A' reduction constants are as follows: C C A, B, C C DJ, DJK, DK, dJ, dK C HJ, HJK, HKJ, HK, hJ, hJK, hK C LJ, LJJK, LJK, LKKJ, LK, lJ, lJK, lKJ, lK C PJ, PJJK, PJK, PKJ, PKKJ, PK, pJ, pJJK, pJK, pKKJ, pK C C The available 'S' reduction constants are as follows: C C A, B, C C DJ, DJK, DK, d1, d2 C HJ, HJK, HKJ, HK, h1, h2, h3 C LJ, LJJK, LJK, LKKJ, LK, l1, l2, l3, l4 C PJ, PJJK, PJK, PKJ, PKKJ, PK, p1, p2, p3, p4, p5 C C C I - number of transitions C N - number of energy levels to be used in fit (if errors are present C in line declarations this may be less than 2*I) C IRED - if set to +-1 then reduction S, otherwise (+-2) reduction A C (positive values representation I.r, negative representation III.r) C K1 - number of fitted constants C FR - information on transitions where C (i,1) observed frequencies C (i,2) calculated frequencies C (i,3) obs-calc C (i,4) relative intensity for blend calculation C ERROR - estimated experimental error in frequency C WEIGHT - weight equal to 1/error**2 C NJA - information on energy levels where C (i,1) J quantum number C (i,2) number of Wang submatrix to which the level belongs =1,2,3,4 for C E+, E-, O+, O- resp. C (i,3) position of the energy level among eigenvalues of the relevant C Wang subm. in ascending order of magnitude C (i,4) number of transition associated with the level, this is negative C for the lower level C IA - pointers to which constants are to be fitted (=1, otherwise 0) C LINFIT - pointers as to what to do with specified transitions C =0 do not fit, predict only C =1 fit C =2 fit intensity weighted mean of several transitions to a single C experimental frequency C JK - quantum numbers for transitions where C (i,1), (i,2), (i,3) are J, K-1, K+1 for upper level resp. C (i,4), (i,5), (i,6) are J, K-1, K+1 for lower level resp. C IPR - default number of iterations C C DESCR - string array holding 8-byte end-of-line annotations C COMBLK - string variable holding all between-the-line comments C concerning transitions C ICOMML(i) - number of transition to be preceded by comment i C ICOMM(i) - position of last character in COMBLK for comment i C ISPACE(i) - number of transition to be preceded by an empty line, the total C number of spacing lines is NSPCES C_____________________________________________________________________________ C C Routines in the order of listing: C C MP C DATAIN - input of data C MODIF - set up internal arrays and online modifications C GENER - generated blocks of lines for specified bands C DATOUT - output of data C DFBET - calculate DFBETA coefficients C hjjset - calculate HJJ coefficients C lzero - adds leading zero to an number in an output string (if necessary) C CONFOR - format constants for printout C ASREDA - set up A reduced Hamiltonian C ASREDS - set up S reduced Hamiltonian C FPLUS - evaluate the Fplus function for the Hamiltonian C FMINUS - evaluate the Fminus function for the Hamiltonian C HDIAG - diagonalise the Hamiltonian matrix C HSORT - sort eigenvalues and eigenvectors out of HDIAG C SORTC - simple ascending order of magnitude sorting of a vector C lentrm - function emulating the LEN_TRIM function C icheck - function to check string for non-FORTRAN characters C C_____________________________________________________________________________ C C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1, lpersc=25) PARAMETER (maxcbl=10000,maxcom=500,maxspa=500) C COMMON /BLKH/H * /BLKR/R * /SORTCC/WORK,IPT COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED * /BLK2/NJA,B * /BLK3/C,D,DD,HJJ * /TBLOCK/COM,T,pinums DIMENSION A(NCONST),B(NCONST),C(NCONST),Z(NCONST),IA(NCONST), * H(MAXDIM,MAXDIM),R(MAXDIM,MAXDIM),DD(NCONST,NCONST) INTEGER*2 NJA(NNJA,4),IPT(NLINES) DIMENSION FR(NLINES,4),JK(NLINES,6),LINFIT(NLINES),WORK(NLINES) REAL*4 D(NLINES,NCONST),EE,ERROR(NLINES),WEIGHT(NLINES), * HJJ(NLINES) character*8 descr(nlines) CHARACTER*6 T(NCONST),TT(NCONST) CHARACTER*1 COM(72) CHARACTER comblk*(maxcbl),pinums(nconst)*3 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c DATA TT/'A =','B =','C =',' DJ =',' DJK =',' DK =', 1 ' dJ =',' dK =','HJ =','HJK =','HKJ =','HK =', 1 'hJ =','hJK =','hK =', 1 ' LJ =',' LJJK=',' LJK =',' LKKJ=',' LK =',' lJ =', 1 ' lJK =',' lKJ =',' lK =', 1 'PJ =','PJJK =','PJK =','PKJ =','PKKJ =','PK =', 1 'pJ =','pJJK =','pJK =','pKKJ =','pK ='/ DATA PINUMS/'100','200','300',' 2',' 11',' 20', 1 '401','410',' 3',' 12',' 21',' 30', 1 '402','411','420', 1 ' 4',' 13',' 22',' 31',' 40','403', 1 '412','421','430', 1 ' 5',' 14',' 23',' 32',' 41',' 50', 1 '404','413','422','431','440'/ DO 1234 L=1,NCONST T(L)=TT(L) 1234 CONTINUE IOUT=0 C IFLAG=0 CALL DATAIN <--- DATAIN 400 CALL MODIF(N,IFLAG) <--- MODIF IF(IOUT.EQ.1)GOTO 41 C I1=0 S=1.D50 IPR=3 C C C START ITERATION C C 14 WRITE(*,6000)SDEV,SWDEV 6000 FORMAT(/' CALCULATION CURRENTLY ON J ---',14X, * 'Last deviations:',F10.4,' (',F6.3,')'/31X,'|'/) I1=I1+1 S1=S DO 15 K=1,I FR(K,2)=0.D0 FR(K,3)=0.D0 DO 15 L=1,NCONST D(K,L)=0.0 15 CONTINUE L1=1 M4=0 WRITE(*,801)NJA(1,1) C C Go through all sorted energy levels C DO 25 K=2,N KK=K-1 IF(NJA(K,1).EQ.NJA(KK,1))GOTO 250 WRITE(*,801)NJA(KK,1) 801 FORMAT(1H+,I4,$) 250 IF(NJA(K,1).EQ.NJA(KK,1).AND.NJA(K,2).EQ.NJA(KK,2))GOTO 25 L2=K-1 N1=NJA(KK,2) J1=NJA(KK,1) C C Set up and diagonalize the appropriate Wang matrix for each group of C common levels C IF(IABS(IRED).EQ.1)THEN CALL ASREDS(M,J1,N1,A,IRED) <--- ASREDS ELSE CALL ASREDA(M,J1,N1,A,IRED) <--- ASREDA ENDIF IEGEN=0 CALL HDIAG(M,IEGEN,K1) <--- HDIAG CALL HSORT(M) <--- HSORT DO 17 L=L1,L2 M=NJA(L,3) M1=NJA(L,4) M2=ABS(M1) FR(M2,2)=FR(M2,2)+SIGN(1,M1)*H(M,M) 17 CONTINUE C C FIND THE DERIVATIVES - since transformation to the space where the Wang C submatrix is diagonal has already been found (eigenvectors), the deri- C vatives for energy levels in this submatrix with respect to a given C constant are obtained by the same similarity transformation of matrix C set up with all other constants set to zero and the required constant C set to 1 (array B) C K1=0 DO 24 J=1,NCONST IF(IA(J).EQ.1)GOTO 18 GOTO 24 18 K1=K1+1 DO 19 L=1,NCONST B(L)=0.D0 19 CONTINUE B(J)=1.D0 IF(IABS(IRED).EQ.1)THEN CALL ASREDS(M3,J1,N1,B,IRED) <--- ASREDS ELSE CALL ASREDA(M3,J1,N1,B,IRED) <--- ASREDA ENDIF DO 124 L=L1,L2 M=NJA(L,3) IF(L.NE.L1.AND.M.EQ.M4)GOTO 23 EE=0.0 DO 22 M1=1,M3 DO 122 M2=M1,M3 N2=2 IF(M1.EQ.M2)N2=1 EE=EE+N2*H(M1,M2)*R(M1,M)*R(M2,M) 122 CONTINUE 22 CONTINUE 23 M1=NJA(L,4) M2=ABS(M1) D(M2,K1)=D(M2,K1)+SIGN(1,M1)*EE M4=M 124 CONTINUE 24 CONTINUE L1=K 25 CONTINUE C C FIND NEW CONSTANTS - inverse of D'D is evaluated by diagonalization using C inv A = U * inv L * U' where U and L contain the eigenvectors and the C eigenvalues of A. C C To improve conditioning D'D is rescaled to D'D.ij/(D'D.ii D'Djj)**1/2 C (scaling factors in Z) C C D'D is placed and diagonalized in the top lh. corner of H, C D'Y is placed in (K1+1)'st column of H C inv(D'D) is placed in array DD C S=0.D0 SW=0.D0 L=0 DO 26 J=1,I C C...unfitted line IF(LINFIT(J).eq.0)then FR(J,3)=DSIGN(FR(J,1),FR(J,2))-FR(J,2) endif C C...single line if(linfit(j).eq.1)then FR(J,3)=DSIGN(FR(J,1),FR(J,2))-FR(J,2) S=S+FR(J,3)**2 SW=SW+(FR(J,3)/ERROR(J))**2 L=L+1 endif C C...blend: deal with the whole blend at the first component if(linfit(j).eq.2)then if(j.eq.1.or.(j.gt.1.and. * ( (fr(J-1,1).ne.fr(j,1)).or.linfit(j-1).ne.2 ) ) )then jfirst=j jlast=j over=fr(j,4)*fr(j,2) below=fr(j,4) do 526 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 527 jlast=jjj over=over+fr(jjj,4)*fr(jjj,2) below=below+fr(jjj,4) 526 continue 527 fr(j,3)=fr(j,1)-over/below do 528 jjj=jfirst+1,i fr(jjj,3)=fr(j,3) if(jjj.eq.jlast)goto 529 528 continue 529 S=S+FR(J,3)**2 SW=SW+(FR(J,3)/ERROR(J))**2 L=L+1 endif endif 26 CONTINUE C NINFIT=L N139=L-K1 IF(N139.LT.0)GOTO 48 IF(N139.GT.0)GOTO 810 WRITE(*,89)L,K1,CHAR(7) SDEV=DSQRT(S) SWDEV=DSQRT(SW) GOTO 811 810 SDEV=DSQRT(S/N139) SWDEV=DSQRT(SW/N139) 811 WRITE(*,'(1x)') WRITE(*,78)I1,SDEV,SWDEV 78 FORMAT(/' ITERATION STEP NO.',I2,11X,'Deviation ', 1 'of fit =',F14.6/24X,'Weighted deviation of fit =',F14.6) C C...in the multiplications below also deal with effect of the whole blend C on encountering the first component C DO 30 J=1,K1 DO 28 K=J,K1 E=0.D0 DO 27 L=1,I IF(LINFIT(L).eq.0)GOTO 27 IF(LINFIT(L).eq.1)E=E+D(L,J)*D(L,K)*WEIGHT(L) IF(LINFIT(L).eq.2)THEN if(L.eq.1.or.(L.gt.1.and. * ( (fr(L-1,1).ne.fr(L,1)).or.linfit(L-1).ne.2 ) ) )then jfirst=L jlast=L over=fr(L,4)*D(L,J)*D(L,K)*WEIGHT(L) below=FR(L,4) do 530 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 531 jlast=jjj over=over+fr(jjj,4)*D(jjj,J)*D(jjj,K)*WEIGHT(jjj) below=below+fr(jjj,4) 530 continue 531 E=E+over/below c E=E+D(L,J)*D(L,K)*WEIGHT(L) endif ENDIF 27 CONTINUE H(J,K)=E 28 CONTINUE Z(J)=SQRT(1.D0/H(J,J)) E=0.D0 DO 29 L=1,I IF(LINFIT(L).eq.0)GOTO 29 IF(LINFIT(L).eq.1)E=E+D(L,J)*FR(L,3)*WEIGHT(L) IF(LINFIT(L).eq.2)THEN if(L.eq.1.or.(L.gt.1.and. * ( (fr(L-1,1).ne.fr(L,1)).or.linfit(L-1).ne.2 ) ) )then jfirst=L jlast=L over=fr(L,4)*D(L,J)*FR(L,3)*WEIGHT(L) below=FR(L,4) do 532 jjj=jfirst+1,i if(linfit(jjj).ne.2.or. * (fr(jjj,1).ne.fr(jfirst,1)) )goto 533 jlast=jjj over=over+fr(jjj,4)*D(jjj,J)*FR(jjj,3)*WEIGHT(jjj) below=below+fr(jjj,4) 532 continue 533 E=E+over/below C E=E+D(L,J)*FR(L,3)*WEIGHT(L) endif ENDIF 29 CONTINUE H(J,K1+1)=E 30 CONTINUE C DO 301 J=1,K1 DO 301 K=J,K1 H(J,K)=H(J,K)*Z(J)*Z(K) H(K,J)=H(J,K) 301 CONTINUE IEGEN=0 C CALL HDIAG(K1,IEGEN,L3) <--- HDIAG C DO 35 J=1,K1 B(J)=0.D0 DO 34 K=1,K1 E=0.D0 DO 33 L=1,K1 E=E+R(J,L)*R(K,L)/H(L,L) 33 CONTINUE E=E*Z(J)*Z(K) DD(J,K)=E B(J)=B(J)+E*H(K,K1+1) 34 CONTINUE C(J)=SQRT(DD(J,J))*SWDEV 35 CONTINUE C WRITE(*,81) 81 FORMAT(/' RESULTS - constant, change, standard deviation :'/) K=0 DO 40 J=1,NCONST IF(IA(J).EQ.1)GOTO 39 GOTO 40 39 K=K+1 A(J)=A(J)+B(K) IF(J.LE.8)THEN WRITE(*,82)T(J),A(J),B(K),C(K) ELSE WRITE(*,1082)T(J),A(J),B(K),C(K) ENDIF 82 FORMAT(1X,A6,3F23.12) 1082 FORMAT(1X,A6,F27.16,2F23.16) 40 CONTINUE C C Convergence tests and output C 403 WRITE(*,'(/'' Do you want to interrupt the fit (1/0) ? '',$)') READ(*,5,ERR=403)IFLAG IF(IFLAG.NE.1)GOTO 6 C 404 WRITE(*,401) 401 FORMAT(1X/ * 9x,'= 0 OBS-CALC + diagnostics = 2 DISK OUTPUT'/ * 9x,'=-1 as above, screen at a time = 3 NEXT ITERATION'/ * 9x,'= 1 MODIFICATIONS = 4 EXIT ASFIT', * 8X,'.... ',$) READ(*,5,ERR=404)IFLAG 5 FORMAT(I5) IF(IFLAG.LT.-1.OR.IFLAG.GT.4)GOTO 404 IF(IFLAG.EQ.1)GOTO 400 IF(IFLAG.EQ.3)GOTO 6 IF(IFLAG.EQ.4)GOTO 99 C c...obs-calc differences c if(iflag.le.0)then lmax=nlines IF(IFLAG.LT.0)lmax=lpersc-1 WRITE(*,'(1X/1x,78(1H-)/1x,t17,''TRANSITION'', * T37,''Frequency'',T52,''Obs-Calc'',T63,''Error'', * T72,''Notes''/1x,78(1H-))') nspace=1 ncomm=2 lcount=3 do 108 j=1,i c if(ispace(nspace).eq.j)then nspace=nspace+1 write(*,'(1x)') lcount=lcount+1 if(lcount.eq.lmax)then write(*,516) read(*,'(i1)',err=515)jj if(jj.eq.1)goto 507 515 lcount=0 endif endif 516 format( * ' Press to continue, 1 to stop .... ',$) c 501 if(icomml(ncomm).eq.j)then write(*,'(1x,a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 lcount=lcount+1 if(lcount.eq.lmax)then write(*,516) read(*,'(i1)',err=506)jj if(jj.eq.1)goto 507 506 lcount=0 goto 501 endif goto 501 endif c IF(LINFIT(J).gt.0)then if(abs(fr(j,3)).le.3.d0*error(j))then if(linfit(j).eq.2) * WRITE(*,93)J,' ',(JK(J,L),L=1,6),' ',FR(J,1),' B', * FR(J,3),ERROR(J),DESCR(J) if(linfit(j).eq.1) * WRITE(*,9)J,(JK(J,L),L=1,6),FR(J,1),FR(J,3),ERROR(J), * DESCR(J) else if(linfit(j).eq.2) * WRITE(*,93)J,'?',(JK(J,L),L=1,6),' ?',FR(J,1),' B', * FR(J,3),ERROR(J),'>3*Error' if(linfit(j).eq.1) * WRITE(*,92)J,(JK(J,L),L=1,6),FR(J,1),FR(J,3),ERROR(J), * '>3*Error' endif else WRITE(*,91)J,(JK(J,L),L=1,6),FR(J,1),FR(J,3), * '--excl' endif lcount=lcount+1 if(lcount.eq.lmax)goto 500 goto 108 c 500 write(*,516) read(*,'(i1)',err=505)jj if(jj.eq.1)goto 507 505 lcount=0 c 108 continue if(iflag.lt.0.and.lcount.gt.lmax-12)then write(*,516) read(*,'(i1)',err=507)jj endif c 507 WRITE(*,78)I1,SDEV,SWDEV c c...most influential lines CALL hjjset(K1,hjjav) <--- hjjset WRITE(*,101) nhj=i-6 if(nhj.lt.1)nhj=1 write(*,102)(ipt(L),hjj(ipt(L)),L=I,nhj,-1) c c...worst fitting lines write(*,105) do 106 L=1,I if(LINFIT(L).gt.0)then work(L)=abs(fr(L,3)/error(L)) if(L.ne.1)then if( (linfit(L).eq.2.and.linfit(L-1).eq.2) .and. * (fr(L,1).eq.fr(L-1,1)) )work(L)=0.d0 endif else WORK(L)=0.d0 endif ipt(L)=L 106 continue int21=1 call sortc(int21,I) <--- sortc nhj=I-6 if(nhj.lt.1)nhj=1 write(*,110)(ipt(L),FR(ipt(L),3)/error(ipt(L)),L=I,nhj,-1) c GOTO 404 ENDIF 9 FORMAT(1X,I4,'. ', 2(I5,2I3), F16.4, F14.4,F8.4,1x,a) 91 FORMAT(1X,I4,'.--',I4,2I3,I5,2I3,'--',F14.4, F14.4,a) 92 FORMAT(1X,I4,'.?', 2(I5,2I3),' ?',F14.4, F14.4,F8.4,1x,a) 93 FORMAT(1X,I4,'.',A, 2(I5,2I3),A, F14.4,a,F12.4,F8.4,1x,a) 101 FORMAT(1x/' MOST INFLUENTIAL TRANSITIONS (line number=hjj):') 102 FORMAT(1x,7(I5,'=0',F4.3)) 105 FORMAT(' WORST FITTING LINES (line number=(o-c)/error):') 110 FORMAT(1x,7(I5,'=',F5.1)) c c...disk output c IF(IFLAG.EQ.2)THEN IOUT=2 CALL DATOUT(SDEV,SWDEV,NINFIT) <--- DATOUT GOTO 404 ENDIF C c...exit to next iteration c 6 IFLAG=1 E=(S1-S)/S IF(DABS(E).LT.0.001D0)GOTO 41 IF(E.LT.-0.1D0)GOTO 41 IF(I1.LT.IPR)GOTO 14 1900 WRITE(*,450) 450 FORMAT(/' MORE CYCLES ? ',$) READ(*,'(I5)',ERR=1900)L IF(L.NE.1)GOTO 41 I1=0 GOTO 14 C 41 CALL DATOUT(SDEV,SWDEV,NINFIT) <--- DATOUT C IF(IFLAG.EQ.1)GOTO 400 GOTO 99 48 WRITE(*,89)L,K1,' - ITERATION STOPPPED' 89 FORMAT(//' ONLY',I3,' MEASURED TRANSITIONS FOR',I3, 1 ' CONSTANTS',A) GOTO 41 C 99 STOP END C C_____________________________________________________________________________ C SUBROUTINE DATAIN IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250) PARAMETER (maxcbl=10000,maxcom=500,maxspa=500) C COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED * /TBLOCK/COM,T,PINUMS DIMENSION JK(NLINES,6),FR(NLINES,4),A(NCONST),IA(NCONST), * LINFIT(NLINES) character*8 descr(nlines) CHARACTER*6 T(NCONST) CHARACTER*1 COM(72) CHARACTER FILNAM*25,cdummy*22,line*90,errmes*50 REAL*4 ERROR(NLINES),WEIGHT(NLINES) CHARACTER comblk*(maxcbl),pinums(nconst)*3 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c IRED=2 C WRITE(*,3344) 3344 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | ASFIT - FITTING OF ASYMMETRIC TOP ROTATIONAL SPECTRUM ', * T79,'|'/ * ' | (reduction A or S, all terms up to decadic)', * T79,'|'/ * ' |',76(1H_),'|'/' version 14.IV.2008',T64,'Zbigniew KISIEL'/) write(*,1910)nlines,jmax 1910 format(' Current limits:',i7,' lines and J up to',i4) 1900 WRITE(*,12) 12 FORMAT(//' DATA TO BE CREATED (+-1), OR READ FROM DISK (+-2,+-3):' 1 //10x,'positive values - reduction A'/ 1 10x,'negative values - reduction S'/ 1 10x,' +-2 - prolate-type representation I.r'/ 1 10x,' +-3 - oblate-type representation III.r'// 1 10x,' ..... ',$) READ(*,14,ERR=1900)IFLAG 14 FORMAT(2I5) C C...IRED=2 -> reduction 'A', IRED=1 -> reduction 'S' C IF(IFLAG.GT.0)GOTO 31 IRED=1 if(iflag.eq.-3)ired=-1 T(7)=' d1 =' T(8)=' d2 =' T(13)='h1 =' T(14)='h2 =' T(15)='h3 =' T(21)='l1 =' T(22)='l2 =' T(23)='l3 =' T(24)='l4 =' T(31)='p1 =' T(32)='p2 =' T(33)='p3 =' T(34)='p4 =' T(35)='p5 =' DO 200 L=21,NCONST A(L)=0.0D0 IA(L)=0.0D0 200 CONTINUE c 31 if(iflag.eq.3)ired=-2 IF(ABS(IFLAG).gt.1)GOTO 15 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C c c...Input from keyboard c 1901 WRITE(*,16) 16 FORMAT(/' COMMENT :') READ(*,17,ERR=1901)COM 17 FORMAT(72A1) 142 write(*,40) 40 format(/' NUMBER OF CONSTANTS TO BE SPECIFIED : ',$) read(*,'(i2)',err=41)nfin if(nfin.lt.1.or.nfin.gt.nconst)goto 142 WRITE(*,18) 18 FORMAT(/' FITTING PARAMETER, INITIAL VALUE OF CONSTANT : ' 1 /) DO 19 L=1,nfin 1903 WRITE(*,20)T(L) 20 FORMAT(1X,A6,2X,$) READ(*,*,ERR=1902)IA(L),A(L) GOTO 19 1902 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 1903 19 CONTINUE I=0 WRITE(*,22) 22 FORMAT(/' LINE ASSIGNMENTS, i.e. six comma separated qn''s, ', 1 'frequency, error: '/ 2 ' (terminate with -1,,,,,,,,)'/) 24 I=I+1 1906 WRITE(*,'('' Line'',i4,'': '',$)')I READ(*,*,ERR=1904)(JK(I,K),K=1,6),FR(I,1),ERROR(I) descr(i)=' ' IF(ERROR(I).EQ.0.)ERROR(I)=1. WEIGHT(I)=1./ERROR(I)**2 GOTO 1905 1904 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 1906 1905 IF(JK(I,1).GE.0)GOTO 24 I=I-1 DO 70 L=1,I LINFIT(L)=1 70 CONTINUE GOTO 41 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C C...input from disk (some features are added so that in DOS compilation there C are no problems with reading UNIX files and files from botched UNIX/DOS C conversion ) C 15 WRITE(*,135) 135 FORMAT(/' NAME OF DATA FILE : ',$) READ(*,'(A)',ERR=15)FILNAM OPEN(2,FILE=FILNAM,ERR=153,STATUS='OLD') c errmes='Bad comment line' READ(2,'(a)',err=152,end=152)line READ(line,3500,ERR=151)COM c errmes='Cannot skip to the number of lines' READ(2,'(a)',err=152,end=152)line 3500 FORMAT(1X,72A1) c errmes='Bad number of transitions' READ(2,'(a)',err=152,end=152)line if(icheck(line).ne.0)goto 157 if(lentrm(line).le.27)line=line(1:lentrm(line))//' ' READ(line(1:27),141,ERR=151)cdummy,I if(I.le.0)then errmes='Zero or negative number of transitions' goto 151 endif c errmes='Bad number of constants' READ(2,'(a)',err=152,end=152)line if(icheck(line).ne.0)goto 157 if(lentrm(line).le.27)line=line(1:lentrm(line))//' ' read(line,141,err=151)cdummy,NCON if(ncon.le.0.or.ncon.gt.nconst)then errmes='Zero, negative or excessive number of constants' goto 151 endif 141 format(a22,i5) errmes='Cannot skip to the list of constants' READ(2,'(a)',err=152,end=152)line c if(ncon.eq.0)ncon=NCONST do 159 k=1,ncon write(errmes,'(a,i3)')'Bad fitting flag or constant number',k READ(2,'(a)',err=152,end=152)line if(icheck(line).ne.0)goto 157 READ(line,150,ERR=151)IA(K),A(K) if(ia(k).lt.0.or.ia(k).gt.1)ia(k)=0 159 continue 150 FORMAT(I5,F30.19) c c...compatibility with files without transitions header c read(2,'(a)',err=151)cdummy if(cdummy(2:11).eq.'----------')then read(2,3500,err=151)cdummy read(2,3500,err=151)cdummy else backspace(2) endif c c...input transitions c nlin=0 nspace=0 ncomm=1 icomm(ncomm)=0 do 3400 L=1,NLINES kkk=nlin+1 c 23 read(2,'(a)',err=241,end=51)line if(icheck(line).ne.0)goto 157 if(lentrm(line).le.70)line=line(1:lentrm(line))//' ' c c...deal with a between-the-line comment if(line(1:1).eq.'!'.or.line(1:1).eq.'%')then if(line(1:1).eq.'%')goto 23 if(ncomm.ge.maxcom)then write(*,35)'Too many comment lines',char(7) goto 23 endif if(icomm(ncomm)+lentrm(line).gt.maxcbl)then write(*,35)'Capacity of the comments buffer exceeded', * char(7) goto 23 endif 35 format(' **** ERROR: ',2a) ncomm=ncomm+1 icomml(ncomm)=L icomm(ncomm)=icomm(ncomm-1)+lentrm(line) comblk(icomm(ncomm-1)+1:icomm(ncomm))=line(1:lentrm(line)) goto 23 endif if(line(1:3).eq.'end')goto 51 if(line(6:10).eq.'-----')goto 23 c c...deal with a separating line if(lentrm(line).lt.20)then nspace=nspace+1 if(nspace.lt.maxspa)ispace(nspace)=L if(nspace.gt.1.and.i.eq.ispace(nspace-1))nspace=nspace-1 goto 23 endif c c...read actual information READ(LINE,140,ERR=241)(JK(kkk,M),M=1,6),FR(kkk,1),ERROR(kkk), * LINFIT(kkk) c read(line(71:80),'(f10.5)',err=3405)fr(kkk,4) goto 3402 3405 fr(kkk,4)=0.d0 c 3402 if(lentrm(line).gt.80)then read(line(81:),'(a)')descr(kkk) else descr(kkk)='' endif c if(linfit(kkk).lt.0.or.linfit(kkk).gt.2)linfit(kkk)=0 140 FORMAT(6I5,F20.6,F15.6,I5) if(linfit(kkk).eq.2.and.fr(kkk,4).eq.0.d0)fr(kkk,4)=1.d0 nlin=nlin+1 3400 CONTINUE 51 I=NLIN if(nlin.eq.nlines)then write(*,3420)nlines 3420 FORMAT(1x/' ***** WARNING: The current maximum of',i7, * ' on the number of lines'/ * ' has been reached.'/ * ' Some lines may not have been read ', * 'and may be lost if'/' the data file is ', * 'overwritten during the fit.'/) PAUSE 'Program can continue but you are advised to rectify this' endif nspces=nspace c CLOSE(2) c c...Check whether blends are specified correctly c nberr=0 nblend=0 bfreq=0.d0 DO 1050 k=1,I c if(linfit(k).ne.2)then if(nblend.eq.0)goto 1050 if(nblend.eq.1)then nberr=nberr+1 if(nberr.eq.1)write(*,'(2x)') write(*,1051)k-1,fr(k-1,1) 1051 format(' **** WARNING: line',i6,F16.5' is only a single', * ' blend component') nblend=0 bfreq=0.d0 goto 1050 endif if(nblend.ge.2)then nblend=0 bfreq=0.d0 goto 1050 endif endif c if(linfit(k).eq.2.and.nblend.eq.0)then nblend=1 bfreq=fr(k,1) goto 1050 endif c if(linfit(k).eq.2.and.nblend.ne.0)then if(fr(k,1).eq.bfreq)then nblend=nblend+1 goto 1050 else if(nblend.eq.1)then write(*,1051)k-1,fr(k-1,1) nberr=nberr+1 nblend=1 bfreq=fr(k,1) goto 1050 else nblend=1 bfreq=fr(k,1) goto 1050 endif endif endif c 1050 continue if(nberr.gt.0)then write(*,1052) 1052 format(1x/' **** WARNING: Program will continue but it is adv', * 'ised to correct the dataset'/15x,'Press ENTER to continue ',$) read(*,'(i5)',err=1053)k 1053 continue endif c c...check for duplicated lines c ndupl=0 DO 1060 kk=1,i-1 if(linfit(kk).eq.0)goto 1060 DO 1062 k=kk+1,i if(linfit(k).eq.0)goto 1062 do 1061 kkk=1,6 if(jk(kk,kkk).ne.jk(k,kkk))goto 1062 1061 continue ndupl=ndupl+1 if(ndupl.eq.1)write(*,'(2x)') write(*,'('' **** WARNING: duplicated lines ''/ * 24x,i5,'':'',3i4,2x,3i4,f16.4/ * 24x,i5,'':'',3i4,2x,3i4,f16.4)') * k,(jk(k,kkk),kkk=1,6),fr(k,1), * kk,(jk(kk,kkk),kkk=1,6),fr(kk,1) 1062 CONTINUE 1060 CONTINUE if(ndupl.gt.0)then write(*,1052) read(*,'(i5)',err=1063)k 1063 continue endif c c...check for split blends c nsplit=0 DO 1055 k=1,I-1 if(fr(k,1).eq.0.d0)goto 1055 if(linfit(k).eq.0)goto 1055 berr=error(k) C DO 1054 kk=k+1,I if(fr(kk,1).eq.0.d0)goto 1054 if(linfit(kk).eq.0)goto 1054 berra=dsqrt(error(kk)**2+berr**2) c if(abs(fr(k,1)-fr(kk,1)).le.berra)then C IF(fr(k,1).eq.fr(kk,1).and.kk-k.eq.1)goto 1054 idiff=0 do 1056 kkk=k+1,kk-1 if(abs(fr(kkk,1)-fr(kk,1)).gt.berra)idiff=idiff+1 1056 continue if(idiff.eq.0)goto 1054 c nsplit=nsplit+1 if(nsplit.eq.1)write(*,'(2x)') write(*,1058) k, (jk(k,j), j=1,6),fr(k,1),error(k), * kk,(jk(kk,j),j=1,6),fr(kk,1),error(kk) endif c 1054 CONTINUE 1055 CONTINUE 1058 format(' **** WARNING - possible split blend: '/ * 9x,i5,':',3i4,2x,3i4,f16.4,f10.4/ * 9x,i5,':',3i4,2x,3i4,f16.4,f10.4) c if(nsplit.gt.0)then write(*,1052) read(*,'(i5)',err=1059)k 1059 continue endif c c...weights c DO 1020 K=1,I IF(ERROR(K).EQ.0.)ERROR(K)=1.0 WEIGHT(K)=1./ERROR(K)**2 1020 CONTINUE GOTO 41 C 153 write(*,'(1x/ '' **** ERROR: nonexistent file, try again'')') goto 15 c 152 write(*,'(1x//'' **** ERROR: '',a// * '' **** END OF FILE encountered on line:''/1x,a//)') * errmes,line(1:lentrm(line)) stop c 151 write(*,'(1x//'' ***** ERROR: '',a// * '' reading line ---''/17x,''|''/1x,a//)') * errmes,line(1:lentrm(line)) if(lentrm(line).eq.0)write(*,'('' (This is an empty line)''//)') stop c 157 write(*,'(1x//'' ***** ERROR: found an embedded non-FORTRAN '', * ''character with ASCII code ='',i4/14x, * ''- please edit this character out of the data set''// * '' reading line ---''/17x,''|''/1x,a//)') * icheck(line),line(1:lentrm(line)) stop c 241 WRITE(*,'(1X/'' ***** ERROR: bad transition number'',I5// * '' reading line ---''/17x,''|''/1x,a//)') * nlin+1,line(1:lentrm(line)) if(lentrm(line).eq.0)write(*,'('' (This is an empty line)''//)') STOP c 243 WRITE(*,'(1X/ * '' ***** ERROR: bad intensity of blend component: '',A// * '' reading line ---''/17x,''|''/1x,a)') * descr(kkk),line STOP C 41 RETURN END C C_____________________________________________________________________________ C SUBROUTINE MODIF(N,IFLAG) C C Modification of data for the fit and identification and sorting of C energy levels: C C N - on output contains number of energy levels identified for fitting C IFLAG - if set to 1 then check of data, listing and reevaluation of NJA C are skipped C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250, * NNJA=2*NLINES+1, MAXDIM=JMAX/2+1, lpersc=25) PARAMETER (maxcbl=10000,maxcom=500,maxspa=500) C COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED * /BLK2/NJA,B * /TBLOCK/COM,T,PINUMS COMMON /BLKH/H * /BLKR/R DIMENSION JK(NLINES,6),FR(NLINES,4),A(NCONST),IA(NCONST), * B(NCONST),LINFIT(NLINES),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM) character*8 descr(nlines) CHARACTER*6 T(NCONST) INTEGER*2 NJA(NNJA,4),IINT(6),IWORK(MAXDIM) CHARACTER*1 COM(72) character wang(4)*2,WANGL*8,pinums(nconst)*3 EQUIVALENCE (WANGL,WANG(1)) REAL*4 ERROR(NLINES),WEIGHT(NLINES) CHARACTER comblk*(maxcbl) dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c NQNS=6 write(wangl,'(a)')'E+E-O+O-' C IF(IFLAG.EQ.1)GOTO 11 C C Check of data, evaluation of NJA (identification of energy levels) and C listing C 29 WRITE(*,30)COM WRITE(*,31) DO 32 K=1,NCONST IF(IA(K).EQ.0.AND.A(K).EQ.0.D0)GOTO 32 WRITE(*,33)T(K),A(K),IA(K) 32 CONTINUE WRITE(*,34)I N=0 NBADL=0 DO 10 K=1,I FR(K,1)=ABS(FR(K,1)) B(1)=JK(K,1)-JK(K,4) C C Test whether J''-J' = 0, -1 or +1 and J less than allowed maximum C IF(ABS(B(1)).GT.1.D0.OR.JK(K,1).GT.JMAX)GOTO 9 DO 6 L=2,5,3 C C Test whether (K-1 + K+1) - J = 0 or 1 C J=JK(K,L)+JK(K,L+1)-JK(K,L-1) IF(J.LT.0.OR.J.GT.1)GOTO 7 N=N+1 NJA(N,1)=JK(K,L-1) NJA(N,4)=K IF(L.EQ.5)NJA(N,4)=-K C C...Wang assignment for representation I.r, old compact coding IF(IRED.GT.0)THEN J=J+(JK(K,L)-JK(K,L)/2*2)*2+1 M=J/2*2 NJA(N,2)=J NJA(N,3)=(JK(K,L)-(J-1)/2-(4-J)/2*M)/2+1 C C...Wang assignment for representation III.r ELSE J=JK(K,L)-JK(K,L)/2*2 IF(JK(K,L-1)/2*2.NE.JK(K,L-1))J=1-J J=J+(JK(K,L+1)-JK(K,L+1)/2*2)*2+1 NJA(N,2)=J IF(JK(K,L-1)-JK(K,L-1)/2*2.EQ.1)THEN IPLUS=J+1 IF(J.GE.3)IPLUS=4-J ELSE IPLUS=J-1 IF(J.GE.3)IPLUS=6-J ENDIF NJA(N,3)=1+(JK(K,L)-JK(K,L+1)+JK(K,L-1)-IPLUS)/4 ENDIF 6 CONTINUE C GOTO 10 C C If error in lower level, eliminate also the upper level from NJA C 7 IF(L.EQ.5)GOTO 8 GOTO 9 8 N=N-1 9 NBADL=NBADL+1 if(NBADL.eq.1)write(*,797) 797 FORMAT(1x/' The following lines have errors in quantum ', * 'numbers and have been'/ * ' automatically excluded from fit:'/) WRITE(*,97)K,(JK(K,M),M=1,NQNS),FR(K,1),ERROR(K) LINFIT(K)=0 10 CONTINUE c 30 FORMAT(/1X,72A1) 31 FORMAT(/' CURRENT CONSTANTS :'/) 33 FORMAT(1X,A6,F28.18,I5) 34 FORMAT(/' TRANSITIONS IN DATA SET:',I6) 97 FORMAT(1X,'Line',I5,': ',2(I5,2I4),F15.4,F10.4, * ' <---INPUT ERROR') c if(nbadl.gt.0)then write(*,1052) 1052 format(1x/' **** WARNING: Program will continue but it is adv', * 'ised to correct the dataset'/15x,'Press ENTER to continue ',$) read(*,'(i5)',err=1053)k 1053 continue endif C C Listing Selection and execution of required modification C 1900 WRITE(*,40) 40 FORMAT(/' SELECT AN OPTION: = -1 list lines page at a time' * /24x, ' = 0 proceed to fit'/ * 24x, ' = 1 list lines'/ * 24x, ' = 2 modifications .... ',$) READ(*,14,ERR=1900)IFLAG 14 FORMAT(I5) IF(IFLAG.EQ.0)GOTO 41 c IF(iabs(IFLAG).eq.1)then lmax=nlines+maxcom+maxspa IF(IFLAG.LT.0)lmax=lpersc-1 WRITE(*,'(1X/1x,78(1H-)/1x,t19,''TRANSITION'', * T39,''Frequency'',T53,''Error'', * T65,''Notes''/1x,78(1H-))') nspace=1 ncomm=2 lcount=3 do 708 j=1,i c if(ispace(nspace).eq.j)then nspace=nspace+1 write(*,'(1x)') lcount=lcount+1 if(lcount.eq.lmax)then write(*,716) read(*,'(i1)',err=711)jj if(jj.eq.1)goto 1900 711 lcount=0 endif endif 716 format( * ' Press to continue, 1 to stop .... ',$) c 701 if(icomml(ncomm).eq.j)then write(*,'(1x,a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 lcount=lcount+1 if(lcount.eq.lmax)then write(*,716) read(*,'(i1)',err=706)jj 706 lcount=0 goto 701 endif goto 701 endif c IF(linfit(J).eq.2) * WRITE(*,790)J,(JK(J,L),L=1,NQNS),FR(J,1),ERROR(J), * '<==Blended line' IF(LINFIT(J).eq.1) * WRITE(*,790)J,(JK(J,L),L=1,NQNS),FR(J,1),ERROR(J),DESCR(J) IF(linfit(J).eq.0) * WRITE(*,790)J,(JK(J,L),L=1,NQNS),FR(J,1),ERROR(J), * '<--Excluded from fit' 790 FORMAT(1X,I4,': ',2(I5,2I4),F15.4,F10.4,2x,A) lcount=lcount+1 if(lcount.eq.lmax)goto 700 goto 708 c 700 write(*,716) read(*,'(i1)',err=705)jj if(jj.eq.1)goto 1900 705 lcount=0 c 708 continue if(iflag.lt.0.and.lcount.gt.lmax-12)then write(*,716) read(*,'(i1)',err=707)jj endif 707 goto 1900 endif c if(Iflag.ne.2)GOTO 1900 C C...M O D I F I C A T I O N M E N U C 11 WRITE(*,70) 70 FORMAT(/' MODIFICATIONS ARE SELECTED AS FOLLOWS :' 1 //' = 0 LISTING AND EXIT TO CALCULATION',T42, 1 '= 6 CHANGE THE COMMENT'/ 1 ' = 1 CHANGE CONSTANTS',T42,'= 7 STORE DATA'/ 1 ' = 2 DELETE LINES',T42,'= 8 EXCLUDE LINES FROM FIT'/ 1 ' = 3 CHANGE LINES',T42,'= 9 INCLUDE LINES IN FIT'/ 1 ' = 4 ADD LINES',T42,'= 10 CHANGE ERROR'/ 1 ' = 5 INSERT LINES',T42,'= 11 GENERATE LINES'/ 1 1X,T42,'= 12 DUMP EIGENVALUES') 72 WRITE(*,71) 71 FORMAT(/' CONTROL PARAMETER (555=EXIT): ',$) READ(*,14,ERR=72)ICONTR if(icontr.eq.555)stop IF(ICONTR.LT.0.OR.ICONTR.GT.12)GOTO 11 C if(icontr.eq. 0)goto 29 exit if(icontr.eq. 1)goto 74 change const if(icontr.eq. 2)goto 75 delete line if(icontr.eq. 3)goto 76 change line if(icontr.eq. 4)goto 77 add line if(icontr.eq. 5)goto 500 insert line if(icontr.eq. 6)goto 178 change cmnt if(icontr.eq. 7)goto 179 store data if(icontr.eq. 8)goto 210 excl line if(icontr.eq. 9)goto 220 incl line if(icontr.eq.10)goto 300 change err if(icontr.eq.11)goto 650 generate if(icontr.eq.12)goto 630 dump c GOTO 72 C C...old: 0 1 2 3 4 5 6 7 8 9 10 11 12 c GOTO(29,74,75,76,77,500,178,179,210,220,300,650,630),ICONTR+1 C C Change of constant (with 1) C 74 WRITE(*,31) LL=NCONST/2 DO 37 L=1,LL WRITE(*,35)(K,T(K),A(K),IA(K),K=L,L+LL,LL) 37 CONTINUE 35 FORMAT(1X,I2,': ',A6,F26.16,I2,2X,I2,': ',A6,F26.16,I2) IF(2*LL.LT.NCONST)WRITE(*,36)NCONST,T(NCONST),A(NCONST), 1 IA(NCONST) 36 FORMAT(39X,2X,I2,': ',A6,F26.16,I2) 1910 WRITE(*,78) 78 FORMAT(/' NUMBER OF CONSTANT TO BE CHANGED (ENTER=exit) ? ',$) READ(*,14,ERR=1910)NN IF(NN.LT.1.OR.NN.GT.NCONST)GOTO 72 1911 WRITE(*,79)T(NN) 79 FORMAT(/' FITTING PARAMETER, NEW VALUE OF ',A6,' ',$) READ(*,*,ERR=1911)IA(NN),A(NN) GOTO 74 C C Line deletion (with 2) C 75 WRITE(*,81) 81 FORMAT(/' NUMBER OF LINE TO BE DELETED ? ',$) READ(*,14,ERR=75)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 I=I-1 DO 1 L=NN,I FR(L,1)=FR(L+1,1) FR(L,4)=FR(L+1,4) LINFIT(L)=LINFIT(L+1) ERROR(L)=ERROR(L+1) WEIGHT(L)=WEIGHT(L+1) DESCR(L)=DESCR(L+1) DO 111 LL=1,NQNS JK(L,LL)=JK(L+1,LL) 111 CONTINUE 1 CONTINUE c do 150 L=1,maxspa if(ispace(L).ge.NN.and.ispace(L).gt.1)ispace(L)=ispace(L)-1 150 continue do 151 L=1,maxcom if(icomml(L).ge.NN.and.icomml(L).gt.1)icomml(L)=icomml(L)-1 151 continue c GOTO 75 C C Change of lines (with 3) C 76 WRITE(*,82) 82 FORMAT(/' NUMBER OF LINE TO BE CHANGED (only frequency if -ve)? ' * ,$) READ(*,14,ERR=76)NN IF(iabs(NN).LT.1.OR.iabs(NN).GT.I)GOTO 72 if(nn.lt.0)then nn=iabs(nn) isele=1 else isele=0 endif 1912 WRITE(*,90)NN,(JK(NN,L),L=1,NQNS),FR(NN,1),ERROR(NN) 90 FORMAT(/ ' LINE NO.',I3,' :',2(I5,2I3),F18.4,F14.4) if(isele.eq.0)then write(*,1913) 1913 format(' NEW LINE : ',$) READ(*,*,ERR=1912)(JK(NN,L),L=1,NQNS),FR(NN,1),ERROR(NN) IF(ERROR(NN).EQ.0.)ERROR(NN)=1. WEIGHT(NN)=1./ERROR(NN)**2 else write(*,1914) 1914 format(23x,' NEW FREQUENCY : ',$) READ(*,*,ERR=1912)FR(NN,1) endif GOTO 76 C C Line addition (with 4) C 77 I=I+1 WRITE(*,190) 190 FORMAT(/' LINE (LINES) TO BE ADDED - qn''s, freq, error:'/ * ' (terminate with -1,,,,,,,,)'/) 191 WRITE(*,'('' Line'',i5,'': '',$)')I READ(*,*,ERR=1915)(JK(I,L),L=1,NQNS),FR(I,1),ERROR(I) GOTO 1916 1915 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 191 1916 IF(JK(I,1).LT.0)GOTO 192 LINFIT(I)=1 FR(I,4)=0.d0 IF(ERROR(I).EQ.0.)ERROR(I)=1. WEIGHT(I)=1./ERROR(I)**2 DESCR(I)=' ' I=I+1 GOTO 191 192 I=I-1 GOTO 72 C C Line insertion (with 5) C 500 WRITE(*,501) 501 FORMAT(/' INSERTIONS ARE TO PRECEDE LINE NO. : ',$) READ(*,14,ERR=500)NINS WRITE(*,502) 502 FORMAT(/' LINE (LINES) TO BE INSERTED - qn''s, freq, error:'/ * ' (terminate with -1,,,,,,,,)'/) 503 WRITE(*,'('' Line'',i5,'': '',$)')NINS READ(*,*,ERR=1918)(IINT(L),L=1,NQNS),FREQ,TEMP GOTO 1919 1918 WRITE(*,'(1H+,A1)')CHAR(7) GOTO 503 1919 IF(IINT(1).lt.0)GOTO 72 I=I+1 DO 504 L=I,NINS,-1 FR(L,1)=FR(L-1,1) FR(L,4)=FR(L-1,4) LINFIT(L)=LINFIT(L-1) ERROR(L)=ERROR(L-1) WEIGHT(L)=WEIGHT(L-1) DESCR(L)=DESCR(L-1) DO 504 LL=1,NQNS JK(L,LL)=JK(L-1,LL) 504 CONTINUE FR(NINS,1)=FREQ FR(NINS,4)=0.d0 DESCR(NINS)=" " LINFIT(NINS)=1 IF(TEMP.EQ.0.)TEMP=1. ERROR(NINS)=TEMP WEIGHT(NINS)=1./TEMP**2 DO 505 LL=1,NQNS JK(NINS,LL)=IINT(LL) 505 CONTINUE NINS=NINS+1 c do 152 L=1,maxspa if(ispace(L).ge.Nins-1)ispace(L)=ispace(L)+1 152 continue do 153 L=1,maxcom if(icomml(L).ge.Nins-1)icomml(L)=icomml(L)+1 153 continue c GOTO 503 C C Change of comment (with 6) C 178 WRITE(*,180)COM 180 FORMAT(/' CURRENT COMMENT :'/1X,A//' NEW COMMENT :') READ(*,181,ERR=178)COM 181 FORMAT(72A1) GOTO 72 C C Saving of data (with 7) C 179 IOUT=1 GOTO 200 C C Line exclusion (with 8) C 210 WRITE(*,211) 211 FORMAT(/' NUMBER OF LINE TO BE EXCLUDED FROM FIT ? ',$) READ(*,14,ERR=210)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 LINFIT(NN)=0 GOTO 210 C C Line inclusion (wiht 9) C 220 WRITE(*,221) 221 FORMAT(/' NUMBER OF LINE TO BE INCLUDED IN FIT ? ',$) READ(*,14,ERR=220)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 LINFIT(NN)=1 GOTO 220 C C Change of error (with 10) C 300 WRITE(*,301) 301 FORMAT(/' ERROR TO BE CHANGED ON LINE NUMBER ? ',$) READ(*,'(i5)',ERR=300)NN IF(NN.LT.1.OR.NN.GT.I)GOTO 72 WRITE(*,302)NN,(JK(NN,K),K=1,NQNS),FR(NN,1),ERROR(NN) 302 FORMAT(/' LINE NO.',I3,' :',2(I5,2I3),F18.4,' +-',F14.4/ 1 ' NEW ERROR : ',$) READ(*,'(f20.10)',ERR=300)TEMP IF(TEMP.EQ.0.)TEMP=1. ERROR(NN)=TEMP WEIGHT(NN)=1./TEMP**2 GOTO 300 C C Generation of lines of a given type for prediction or as templates C for input (with 11) C 650 write(*,651) 651 format(/' SELECT 0 = generation according to band rules ', * '(and exit)'/ * ' 1 = generation by cloning a given line'// * ' ..... ',$) READ(*,14,ERR=650)ICLONE if(iclone.eq.0)goto 600 if(iclone.ne.1)goto 650 c c...clones of a supplied line c iold=i CALL GENER(ICLONE,IBAND,IBOTH) goto 652 c C...bands c 600 WRITE(*,601) 601 FORMAT(/' SELECT TYPE OF BAND (-ve parameter specifies only ', * 'one of two like components):'// * T10,'0 = EXIT this option'/ * T10,'1 = type- I(-) a-R.0,1 (various K-1 for given J)'/ * T9,' 2 = type- I(+) c-R.1,0 (various K+1 for given J)'/ * T10,'3 = type-II(+) a-R.0,1 \ J decreases down ', * 'the band,',/ * T10,'4 = type-II(+) b-R.1,1 + b.R.-1,1 / (K+1 large and', * ' degenerate)'/ * T10,'5 = type-II(-) b-R.1,1 + b.R.-1,1 \ J increases down ', * 'the band,'/ * T10,'6 = type-II(-) c-R.1,0 / (K-1 large and', * ' degenerate)'/ * T10,'7 = a-Q.1,0'/ * T10,'8 = b-Q.1,-1 + b-Q.-1,1'/ * T10,'9 = c-Q.1,0 + c-Q.1,-2'/ * T9,'10 = c-R.-1,2 '/ * T40,'.... ',$) IBOTH=1 READ(*,14,ERR=600)IBAND IF(IBAND.EQ.0)GOTO 72 IF(IBAND.LT.0)THEN IBOTH=0 IBAND=IABS(IBAND) ENDIF IF(IBAND.LT.1.OR.IBAND.GT.11)GOTO 600 IOLD=I C CALL GENER(ICLONE,IBAND,IBOTH) C 652 WRITE(*,606)I-IOLD,I 606 FORMAT(1X/I5,' lines added, total number of lines is now:',I5) GOTO 650 C C...Dump of energy levels for a specified J+1<-J transition (with 12) C C NOTE: K-1,K+1 parity: C C J.even J.odd J.even J.odd C E+ I: ee eo III: ee oe C E- eo ee oe ee C O+ oo oe eo oo C O- oe oo oo eo C C K-1 = ( Tau + J+1 ) /2 C K+1 = ( J+1 - Tau ) /2 C C (from Table 7.5 of Gordy and Cook, p.243 third ed. these parities are C for I.l and III.l and in that table entries for I.r and III.r differ C by reversal of O+ and O- lines. --To be resolved!) C 630 NEWDMP=1 636 WRITE(*,631) 631 FORMAT(1X/' J for J+1<-J transition: ',$) READ(*,14,ERR=636)J IF(J.EQ.0.or.J.eq.555.or.j-1.gt.jmax)GOTO 72 OPEN(1,FILE='DUMP.ASF',ACCESS='APPEND') IF(NEWDMP.EQ.1)THEN WRITE(1,30)COM WRITE(1,31) LL=NCONST/2 DO 635 L=1,LL WRITE(1,35)(K,T(K),A(K),IA(K),K=L,L+LL,LL) 635 CONTINUE IF(2*LL.LT.NCONST)WRITE(1,36)NCONST,T(NCONST),A(NCONST), * IA(NCONST) ENDIF C IEGEN=1 write(*,634)j 634 format(' ---- Now working on J=',i4) do 640 nwang=1,4 if(IABS(ired).ne.1)then call asreda(l,j,nwang,a,IRED) else call asreds(l,j,nwang,a,IRED) endif CALL HDIAG(L,IEGEN,NIT) CALL HSORT(L) WRITE(1,632)J,wang(nwang) L=0 IF(IRED.GT.0)THEN IPLUS=5-NWANG IF(NWANG.EQ.1)IPLUS=0 ELSE IF(J-J/2*2.EQ.1)THEN IPLUS=NWANG+1 IF(NWANG.GE.3)IPLUS=4-NWANG ELSE IPLUS=NWANG-1 IF(NWANG.GE.3)IPLUS=6-NWANG ENDIF ENDIF DO 642 LL=-J+IPLUS,J,4 L=L+1 IWORK(L)=LL 642 CONTINUE WRITE(1,633)(J,(J+1+IWORK(LL))/2,(J+1-IWORK(LL))/2, * H(LL,LL),LL=1,L) 640 continue 632 FORMAT(1X/' J=',I3,', Wang matrix ',A,':') 633 FORMAT(3(I6,',',I2,',',I2,F14.3)) C J=J+1 write(*,634)j do 641 nwang=1,4 if(IABS(ired).ne.1)then call asreda(l,j,nwang,a,IRED) else call asreds(l,j,nwang,a,IRED) endif CALL HDIAG(L,IEGEN,NIT) CALL HSORT(L) WRITE(1,632)J,wang(nwang) L=0 IF(IRED.GT.0)THEN IPLUS=5-NWANG IF(NWANG.EQ.1)IPLUS=0 ELSE IF(J-J/2*2.EQ.1)THEN IPLUS=NWANG+1 IF(NWANG.GE.3)IPLUS=4-NWANG ELSE IPLUS=NWANG-1 IF(NWANG.GE.3)IPLUS=6-NWANG ENDIF ENDIF DO 643 LL=-J+IPLUS,J,4 L=L+1 IWORK(L)=LL 643 CONTINUE WRITE(1,633)(J,(J+1+IWORK(LL))/2,(J+1-IWORK(LL))/2, * H(LL,LL),LL=1,L) 641 continue CLOSE(1) NEWDMP=0 GOTO 636 C C------------------------------------- C C Sort of energy levels in NJA according to J, Wang submatrix and order C of eigenvalue in the submatrix C 41 DO 13 J=1,N-1 K1=J+1 L1=J DO 12 K=K1,N IF(NJA(K,1)-NJA(L1,1))411,401,12 401 IF(NJA(K,2)-NJA(L1,2))411,402,12 402 IF(NJA(K,3)-NJA(L1,3))411,12,12 411 L1=K 12 CONTINUE IF(L1.EQ.J)GOTO 13 DO 400 L=1,4 M=NJA(J,L) NJA(J,L)=NJA(L1,L) NJA(L1,L)=M 400 CONTINUE 13 CONTINUE N=N+1 NJA(N,1)=0 NJA(N,2)=0 NJA(N,3)=0 C 200 RETURN END C C_____________________________________________________________________________ C SUBROUTINE GENER(ICLONE,IBAND,IBOTH) C C GENERATION OF TRANSITIONS IN BANDS WHERE C C ICLONE =1 lines are to be cloned from supplied quantum numbers, in which C case IBAND and IBOTH do not apply C =0 lines are to be generated from band rules C IBAND = band type as below or in header in MODIF c IBOTH=1 both transitions which may become degenerate for each line to be C calculated C IBOTH=0 only one of the like transitions to be calculated C C Generated transitions will be assigned the error carried by the last C transition in the data set C C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250) PARAMETER (maxcbl=10000,maxcom=500,maxspa=500) C COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED COMMON /annot/ispace,icomm,icomml,comblk,nspces character*8 descr(nlines) REAL*4 ERROR(NLINES),WEIGHT(NLINES) DIMENSION JK(NLINES,6),FR(NLINES,4),A(NCONST),IA(NCONST), * LINFIT(NLINES) CHARACTER comblk*(maxcbl) dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) C C C___LINE CLONING C iold=i if(iclone.eq.0)goto 700 702 write(*,701) 701 format(1x/' Type in six quantum numbers for the clone template'/ * ' (comma separated)'// * ' ..... ',$) read(*,*,err=702)ju,kau,kcu,jl,kal,kcl if(iabs(ju-jl).gt.1)goto 702 n=iabs(ju-(kau+kcu)) if(n.lt.0.or.n.gt.1)goto 702 n=iabs(jl-(kal+kcl)) if(n.lt.0.or.n.gt.1)goto 702 kaddu=kau+kcu-ju kaddl=kal+kcl-jl c 704 write(*,703) 703 format(1x/' Select 1 = fix J'/ * ' 2 = fix K-1'/ * ' 3 = fix K+1'// * ' (-ve value selects both possible', * ' components)'// * ' ..... ',$) read(*,*,err=704)nfix iboth=0 if(nfix.lt.0)then nfix=-nfix iboth=1 endif if(nfix.lt.1.or.nfix.gt.3)goto 704 c 706 write(*,705) 705 format(1x/' Shift in quantum number relative to this transition'/ * ' (comma separated)'// * ' ..... ',$) read(*,*,err=706)minqn,maxqn if(minqn.gt.maxqn)goto 706 c if(iboth.eq.1)then if(nfix.eq.1)then 716 write(*,715) 715 format(1x/' Select -1 = prolate (K-1) doubling'/ * ' 1 = oblate (K+1) doubling'// * ' ..... ',$) read(*,*,err=716)kdubl if(iabs(kdubl).ne.1)goto 716 else if(nfix.eq.2)kdubl=-1 if(nfix.eq.3)kdubl=1 endif endif c C...Generation loop C do 707 k=minqn,maxqn c c...fixed J c if(nfix.eq.1)then jjl=jl kkl=kal+k if(kkl.gt.jjl)goto 707 kkkl=jjl-kkl+kaddl if(kkkl.gt.jjl)goto 707 jju=ju kku=kau+k if(kku.gt.jju)goto 707 kkku=jju-kku+kaddu if(kkku.gt.jju)goto 707 endif c c...fixed K-1 c if(nfix.eq.2)then jju=ju+k if(jju.lt.0)goto 707 jjl=jl+k if(jjl.lt.0)goto 707 kku=kau kkl=kal kkku=jju-kku+kaddu kkkl=jjl-kkl+kaddl endif c c...fixed K+1 c if(nfix.eq.3)then jju=ju+k if(jju.lt.0)goto 707 jjl=jl+k if(jjl.lt.0)goto 707 kkku=kcu kkkl=kcl kku=jju-kkku+kaddu kkl=jjl-kkkl+kaddl endif c if(kkl .lt.0.or.kku .lt.0)goto 707 if(kkkl.lt.0.or.kkku.lt.0)goto 707 i=i+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(I-1) WEIGHT(I)=1.D0/ERROR(I)**2 DESCR(I)=' ' LINFIT(I)=0 c JK(I,1)=jju JK(I,4)=jjl JK(I,2)=kku JK(I,5)=kkl JK(I,3)=kkku JK(I,6)=kkkl c if(iboth.eq.1)then if(kdubl.eq.-1.and.( kkl.eq.0.or. kku.eq.0))goto 707 if(kdubl.eq. 1.and.(kkkl.eq.0.or.kkku.eq.0))goto 707 c if(kdubl.eq.-1)then if(kku+kkku.eq.jju)then kkku=kkku+1 else kkku=kkku-1 endif if(kkl+kkkl.eq.jjl)then kkkl=kkkl+1 else kkkl=kkkl-1 endif else if(kku+kkku.eq.jju)then kku=kku+1 else kku=kku-1 endif if(kkl+kkkl.eq.jjl)then kkl=kkl+1 else kkl=kkl-1 endif endif c i=i+1 FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(I-1) WEIGHT(I)=1.D0/ERROR(I)**2 DESCR(I)=' ' LINFIT(I)=0 JK(I,1)=jju JK(I,4)=jjl JK(I,2)=kku JK(I,5)=kkl JK(I,3)=kkku JK(I,6)=kkkl endif c 707 continue return C C C___GENERATION ACCORDING TO BAND RULES C C...a-R.0,1, type-I(-): C 700 IERROR=I IF(IBAND.EQ.1)THEN 603 WRITE(*,610)'Type-I(-) a-R.0,1: J' 610 FORMAT(1X/1X,8(1H-),T10,A,' = ',$) 602 FORMAT(1X/T10,A,' = ',$) READ(*,14,ERR=603)JBAND 14 FORMAT(I5) IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 603 604 WRITE(*,602)' K-1 range (lower,upper)' READ(*,*,ERR=604)KMLOW,KMHIGH IF(KMHIGH.LT.KMLOW)GOTO 604 IF(KMLOW.LT.0.OR.KMLOW.GT.JBAND)GOTO 604 IF(KMHIGH.GT.JBAND)GOTO 604 DO 605 K=KMLOW,KMHIGH DO 605 KK=0,1 IF((K.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 605 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=JBAND+1 JK(I,4)=JBAND JK(I,2)=K JK(I,5)=K JK(I,3)=JK(I,1)-JK(I,2)+KK JK(I,6)=JK(I,4)-JK(I,5)+KK FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 DESCR(I)=' ' LINFIT(I)=0 605 CONTINUE ENDIF C C...c-R.1,0, type-I(+): C IF(IBAND.EQ.2)THEN 653 WRITE(*,610)'Type-I(+) c-R.1,0: J range (lower,upper)' READ(*,*,ERR=653)JBLOW,JBHIGH IF(JBHIGH.LT.JBLOW)GOTO 653 IF(JBLOW.LT.0)GOTO 653 IF(JBHIGH.GT.JMAX)GOTO 653 654 WRITE(*,602)' K+1 range (lower,upper)' READ(*,*,ERR=654)KPLOW,KPHIGH IF(KPHIGH.LT.KPLOW)GOTO 654 IF(KPLOW.LT.0.OR.KPLOW.GT.JBHIGH)GOTO 654 IF(KPHIGH.GT.JBHIGH)GOTO 654 DO 655 JBAND=JBLOW,JBHIGH DO 655 K=KPLOW,KPHIGH IF(K.GT.JBAND)GOTO 655 DO 656 KK=0,1 IF((K.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 656 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=JBAND+1 JK(I,4)=JBAND JK(I,3)=K JK(I,6)=K JK(I,2)=JK(I,1)-JK(I,3)+KK JK(I,5)=JK(I,4)-JK(I,6)+KK FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 LINFIT(I)=0 DESCR(I)=' ' 656 CONTINUE 655 CONTINUE ENDIF C C...a-R.0,1 and b-R, type-II(+): bands are formed when (A+B)=(l/m)(2C), where C l and m are integers such that (l/m)>1 C IF(IBAND.EQ.3.OR.IBAND.EQ.4)THEN 616 WRITE(*,610) 'Type-II(+): l and m defining n=(l/m)' read(*,*,err=616)LBAND,MBAND if(lband.lt.2.or.mband.lt.1)goto 616 if(lband.le.mband)goto 616 IF(IBAND.EQ.3)THEN 611 WRITE(*,610)'Type-II(+) a-R.0,1: J of leading line' READ(*,14,ERR=611)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 611 ELSE 615 WRITE(*,610)'Type-II(+) b-R: J of leading line' READ(*,14,ERR=615)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 615 ENDIF isband=0 if(mband.gt.1)then write(*,'(1x/'' Subbands are possible for following K-1: '', * 15i3)')(kk,kk=0,mband-1) 630 write(*,602)' K-1 for leading line in subband' read(*,14,err=630)isband if(isband.lt.0.or.isband.gt.mband-1)goto 630 endif c 612 WRITE(*,602) ' K-1 range (lower,upper)' READ(*,*,ERR=612)KMLOW,KMHIGH IF(KMHIGH.LT.KMLOW)GOTO 612 IF(KMLOW.LT.0.OR.KMLOW.GT.JBAND)GOTO 612 IF(KMHIGH.GT.JBAND)GOTO 612 jline=jband-(mband-lband) DO 6614 K=ISBAND,KMHIGH,MBAND jline=jline+(mband-lband) if(k.lt.kmlow)goto 6614 if(jline-k.lt.0)goto 6614 DO 614 KK=0,1 IF(IBOTH.EQ.0.AND.KK.EQ.1)GOTO 6614 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,4)=jline JK(I,1)=JK(I,4)+1 IF(IBAND.EQ.3)THEN aR 0,1 JK(I,2)=K+KK JK(I,5)=K+KK JK(I,6)=JLINE-K JK(I,3)=JK(I,6)+1 ELSE bR 1,1 JK(I,5)=K+KK bR 1,-1 JK(I,2)=K+1-KK JK(I,6)=JLINE-K JK(I,3)=JK(I,6)+1 ENDIF FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 LINFIT(I)=0 DESCR(I)=' ' 614 CONTINUE 6614 CONTINUE ENDIF C C...b-R and c.R.1,0, type-II(-): bands are formed when 2A=(l/m)(B+C), where C l and m are integers such that (l/m)>1 C IF(IBAND.EQ.5.OR.IBAND.EQ.6)THEN 516 WRITE(*,610) 'Type-II(-): l and m defining n=(l/m)' read(*,*,err=516)LBAND,MBAND if(lband.lt.2.or.mband.lt.1)goto 516 if(lband.le.mband)goto 516 c IF(IBAND.EQ.5)THEN 511 WRITE(*,610)'Type-II(-) b-R: J of leading line' READ(*,14,ERR=511)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 511 ELSE 515 WRITE(*,610)'Type-II(-) c-R.1,0: J of leading line' READ(*,14,ERR=515)JBAND IF(JBAND.LT.0.OR.JBAND.GT.JMAX)GOTO 515 ENDIF c write(*,'(1x/'' Subbands are possible for following K+1: '', * 15i3)')(kk,kk=0,lband-1) 530 write(*,602)' K-1 for leading line in subband' read(*,14,err=530)isband if(isband.lt.0.or.isband.gt.lband-1)goto 530 c 512 WRITE(*,602) ' K+1 range (lower,upper)' READ(*,*,ERR=512)KMLOW,KMHIGH IF(KMHIGH.LT.KMLOW)GOTO 512 IF(KMLOW.LT.0.OR.KMLOW.GT.JBAND)GOTO 512 IF(KMHIGH.GT.JBAND)GOTO 512 jline=jband-(lband-mband) DO 5514 K=isband,KMHIGH,lband jline=jline+(lband-mband) if(k.lt.kmlow)goto 5514 if(jline-k.lt.0)goto 5514 DO 514 KK=0,1 IF(IBOTH.EQ.0.AND.KK.EQ.1)GOTO 514 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,4)=jline JK(I,1)=JK(I,4)+1 IF(IBAND.EQ.6)THEN cR 1,0 JK(I,6)=K+KK JK(I,3)=K+KK JK(I,5)=JLINE-K JK(I,2)=JK(I,5)+1 ELSE bR 1,1 JK(I,6)=K+KK bR 1,-1 JK(I,3)=K+1-KK JK(I,5)=JLINE-K JK(I,2)=JK(I,5)+1 ENDIF FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 LINFIT(I)=0 DESCR(I)=' ' 514 CONTINUE 5514 CONTINUE ENDIF C C...a-Q.1,0 C IF(IBAND.EQ.7)THEN WRITE(*,'(1X/'' ***** SORRY, NOT YET IMPLEMENTED''/)') ENDIF C C...b-Q.1,-1 + b-Q.-1,1 C IF(IBAND.EQ.8)THEN 607 WRITE(*,610)'b-Q.1,-1 (b-Q.-1,1): K-1' READ(*,14,ERR=607)KBAND IF(KBAND.LT.0.OR.KBAND.GT.JMAX)GOTO 607 608 WRITE(*,602)' J range (lower,upper)' READ(*,*,ERR=608)JLOW,JHIGH IF(JHIGH.LT.JLOW)GOTO 608 IF(JLOW.LT.KBAND)GOTO 608 IF(JHIGH.GT.JMAX)GOTO 608 DO 609 J=JLOW,JHIGH DO 609 KK=0,1 IF((KBAND.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 609 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=J JK(I,4)=J JK(I,2)=KBAND+1 JK(I,5)=KBAND JK(I,6)=JK(I,4)-JK(I,5)+KK JK(I,3)=JK(I,6)-1 FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 LINFIT(I)=0 DESCR(I)=' ' 609 CONTINUE ENDIF C C...c-Q.1,0 + c-Q.1,-2 C IF(IBAND.EQ.9)THEN 670 WRITE(*,610)'c-Q.1,0 (c-Q.1,-2): K-1' READ(*,14,ERR=670)KBAND IF(KBAND.LT.0.OR.KBAND.GT.JMAX)GOTO 670 671 WRITE(*,602)' J range (lower,upper)' READ(*,*,ERR=671)JLOW,JHIGH IF(JHIGH.LT.JLOW)GOTO 671 IF(JLOW.LT.KBAND.OR.JLOW.EQ.0)GOTO 671 IF(JHIGH.GT.JMAX)GOTO 671 DO 672 J=JLOW,JHIGH DO 672 KK=0,1 IF((KBAND.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 672 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=J JK(I,4)=J JK(I,2)=KBAND+1 JK(I,5)=KBAND IF(KK.EQ.0)THEN JK(I,6)=JK(I,4)-JK(I,5) JK(I,3)=JK(I,6) ELSE JK(I,6)=JK(I,4)-JK(I,5)+1 JK(I,3)=JK(I,6)-2 ENDIF FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 LINFIT(I)=0 DESCR(I)=' ' IF(JK(I,2)+JK(I,3).GT.JK(I,1)+1)I=I-1 IF(JK(I,2).LT.0)I=I-1 672 CONTINUE ENDIF C C...c-R.-1,2 C IF(IBAND.EQ.10)THEN 650 WRITE(*,610)'c-R.-1,2: K-1' READ(*,14,ERR=650)KBAND IF(KBAND.LT.1.OR.KBAND.GT.JMAX-2)GOTO 650 651 WRITE(*,602)' J range (lower,upper)' READ(*,*,ERR=651)JLOW,JHIGH IF(JHIGH.LT.JLOW)GOTO 651 IF(JLOW.LT.KBAND)GOTO 651 IF(JHIGH.GT.JMAX)GOTO 651 DO 652 J=JLOW,JHIGH DO 652 KK=0,1 IF((KBAND.EQ.0.OR.IBOTH.EQ.0).AND.KK.EQ.1)GOTO 652 I=I+1 if(i-1.eq.iold)then nspces=nspces+1 ispace(nspces)=i endif JK(I,1)=J+1 JK(I,4)=J JK(I,2)=KBAND+KK-1 JK(I,5)=KBAND+KK JK(I,6)=JK(I,4)-JK(I,5)+KK JK(I,3)=JK(I,6)+2 FR(I,1)=0.0D0 FR(I,4)=0.0D0 ERROR(I)=ERROR(IERROR) WEIGHT(I)=1.D0/ERROR(I)**2 LINFIT(I)=0 DESCR(I)=' ' IF(JK(I,2)+JK(I,3).GT.JK(I,1)+1)I=I-1 IF(JK(I,2).LT.0)I=I-1 652 CONTINUE ENDIF C RETURN END C C_____________________________________________________________________________ C SUBROUTINE DATOUT(SDEV,SWDEV,NINFIT) C C OUTPUT OF DATA AND RESULTS FROM ASFIT - BOTH TO SCREEN AND TO DISK C C SDEV - standard deviation of the fit C SWDEV - deviation per unit weight C NINFIT - number of different frequency transitions in fit C NCONST - total no of fittable constants C also: C IOUT - if set to 1 on input, then only disk output of data is carried out C if set to 2, then output menu is provided for selection of C disk output of data and/or intermediate results C K - output mode which is determined from the menu at the bottom of C the routine C =0 no output C =1 standard output (including c.c., worst and most infl.transitions) C =2 above + c.c. distribution plot C =3 above + contributions or hjj/DFBETAS C =4 constants file for ASROT C =5 frequency listing in deposition format C =6 .LIN file for SPFIT C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250, MAXDIM=JMAX/2+1, * lpersc=25) PARAMETER (maxcbl=10000,maxcom=500,maxspa=500) C COMMON /BLKH/H * /BLKR/R * /SORTCC/WORK,IPT COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED * /BLK3/C,D,DD,HJJ * /TBLOCK/COM,T,PINUMS INTEGER*2 JK(NLINES,6),IA(NCONST),TPT(NCONST),LINFIT(NLINES), * UNITS(NCONST),NCC(20),namcd(10),tptinv(nconst),IPT(NLINES) REAL*8 FR(NLINES,4),A(NCONST),C(NCONST),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM),DD(NCONST,NCONST),WORK(NLINES), * DFBETA(NCONST),PICON(NCONST) REAL*4 D(NLINES,NCONST),DER,DIV,ERROR(NLINES),WEIGHT(NLINES), * HJJ(NLINES) character*8 descr(nlines),tdescr CHARACTER T(NCONST)*6,linout*500,lincom*100,linrel*10,annot*10 CHARACTER*25 FILNAM,conval,erval CHARACTER*38 OUTCON(NCONST) CHARACTER*7 UNIT,NAMUN(6) CHARACTER COM(72),FORMAR(54),NAMRED,namrep*5,cdum1*20,cdum2*20 CHARACTER*54 FORMAL,F1,F2 EQUIVALENCE (FORMAL,FORMAR(1)) CHARACTER comblk*(maxcbl),csymb,pinums(nconst)*3,vv*2 dimension ispace(maxspa),icomm(maxcom),icomml(maxcom) COMMON /annot/ispace,icomm,icomml,comblk,nspces c iexdep=0 NOUT=0 C C...UNITS defines the units in which various constants are to be printed C out, NAMUN contains the names of these units DATA NAMUN/'MHz ','kHz ',' Hz ','mHz ','microHz', * 'nHz '/ DATA UNITS/1,1,1, 2,2,2,2,2, 3,3,3,3,3,3,3, 4,4,4,4,4,4,4,4,4, * 5,5,5,5,5,5,5,5,5,5,5/ C C...F1 - format for first line of printout of contributions, F2 - format C for all following lines for the same transition F1='( I4,3H. ,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' F2='(1X,6H ,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' namrep='I.r ' NAMRED='A' IF(iabs(IRED).EQ.1)NAMRED='S' IF(IRED.LT.0)NAMREP='III.r' IF(IOUT.EQ.1)GOTO 100 C K1=0 DO 9 NN=1,NCONST IF(IA(NN).NE.1)GOTO 9 K1=K1+1 TPT(K1)=NN TPTINV(nn)=k1 9 CONTINUE C C...go to the output menu C IF(IOUT.EQ.2)GOTO 300 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...Output of results (standard listing) which can be routed to screen C (NOUT=0) or to disk (NOUT=3) C 25 if(k.ne.5)then WRITE(NOUT,1)COM,NAMRED,NAMREP else write(nout,61)com,namred,namrep endif 1 FORMAT( T54,'ASFIT version 14.IV.2008'/ * 79(1H_)//1x,72A1/79(1H_)//' Fit of Watson''s ', * 'Hamiltonian --> reduction-',A,', representation-',A,':'//) 61 FORMAT(79(1H_)//1X,72A1/79(1H_)//' Obs-calc differences', * ' are from fit of Watson''s reduced Hamiltonian in'/33x, * '---> reduction-',A,', representation-',A,'<---'/) L=0 IF(K.NE.5)THEN DO 2 N=1,NCONST IF(A(N).EQ.0.D0)GOTO 2 IF(IA(N).EQ.0)GOTO 3 L=L+1 WRITE(NOUT,4)T(N),A(N),C(L) GOTO 2 3 WRITE(NOUT,5)T(N),A(N) 2 CONTINUE 4 FORMAT(1X,A6,F30.19,' +-',F25.19) 5 FORMAT(1X,A6,F30.19) WRITE(NOUT,405)SDEV,SWDEV 405 FORMAT(/' Deviation of fit = ',F14.6/ 1 ' Weighted deviation of fit = ',F14.6) WRITE(NOUT,'(1x)') if(nout.eq.0)then 902 write(*,901) 901 format(' F I T C O M P L E T E D: ', 1 '0 = continue'/ 1 29x,' 1 = scroll through fitted lines'/ 1 29x,'-1 = as above, screen at a time .... ',$) read(*,'(I2)',err=902)nscrol if(iabs(nscrol).ne.1)goto 300 WRITE(NOUT,'(1x)') WRITE(NOUT,53) WRITE(NOUT,6) 6 FORMAT(17X,'TRANSITION',T42,'Observed',T56,'Obs-Calc',T68, 1 'Error') WRITE(NOUT,53) else WRITE(NOUT,553) WRITE(NOUT,666) 666 FORMAT(17X,'TRANSITION',T42,'Observed',T56,'Obs-Calc',T68, 1 'Error Note Blend components') WRITE(NOUT,553) endif ELSE WRITE(NOUT,67) 67 FORMAT(/10X,'TRANSITION',T35,'Observed',T47,'Obs-Calc',T59, 1 'Error Note wt for'/T76,'blend'/) ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...TRANSITIONS, listed for all k but: C C if k=5 then deposition type printout C if NOUT=0 and nscrol=-1 then scrolled page at a time to screen c nspace=1 ncomm=2 if(iabs(nscrol).eq.1)then lmax=nlines+maxcom+maxspa IF(nscrol.LT.0)lmax=lpersc-1 lcount=3 endif DO 7 N=1,I lendes=lentrm(descr(n)) if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') lcount=lcount+1 if(nout.eq.0.and.lcount.eq.lmax)then write(*,716) read(*,'(i1)',err=707)jj if(jj.eq.1)goto 709 707 lcount=0 endif endif 716 format( * ' Press to continue, 1 to stop .... ',$) c 33 if(icomml(ncomm).eq.n)then if(nout.eq.0)then write(NOUT,'(1x,a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) else write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) endif ncomm=ncomm+1 lcount=lcount+1 if(nout.eq.0.and.lcount.eq.lmax)then write(*,716) read(*,'(i1)',err=706)jj if(jj.eq.1)goto 709 706 lcount=0 goto 33 endif goto 33 endif c c...transition in fit c if(lentrm(descr(n)).eq.1.and.descr(n)(1:1).eq.' ')then annot=' ' else annot=' /'//descr(n) endif c IF(LINFIT(N).gt.0)then c c...non-deposition if(k.ne.5)then if(abs(fr(n,3)).le.3.d0*error(n))then if(nout.ne.0)then if(linfit(n).eq.1) * WRITE(NOUT,81)N,' ',(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERROR(N),annot if(linfit(n).eq.2) * WRITE(NOUT,181)N,'B',(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERROR(N),annot,fr(n,4),FR(n,2) else WRITE(NOUT,721)N,' ',(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERROR(N) endif else csymb='?' if(nout.ne.0)then if(linfit(n).eq.1) * WRITE(NOUT,81)N,csymb,(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERROR(N),' >3*Err ',annot if(linfit(n).eq.2) * WRITE(NOUT,181)N,csymb,(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERROR(N),' >3*Err ',fr(n,4),FR(n,2) else tdescr='>3*Err' if(abs(fr(n,3)).le.3.d0*error(n))tdescr='<=Blnd' WRITE(NOUT,721)N,csymb,(JK(N,L),L=1,6),FR(N,1),FR(N,3), * ERROR(N),tdescr(1:6) endif endif c c...deposition else ndig=10 if(error(n).gt.0.01d0)then write(cdum1,'(f10.3)')fr(n,3) write(cdum2,'(f10.3)')error(n) else write(cdum1,'(f10.4)')fr(n,3) write(cdum2,'(f10.4)')error(n) endif call lzero(cdum1,ndig) call lzero(cdum2,ndig) if(fr(n,4).eq.0.d0)then if(error(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),annot else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),annot endif else if(error(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),annot,fr(n,4) else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),annot,fr(n,4) endif endif endif c c...transition excluded from fit c c...non-deposition else if(k.ne.5)then if(nout.ne.0)then WRITE(NOUT,81)N,'E',(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERROR(n),'--excl ',annot else WRITE(NOUT,721)N,'E',(JK(N,L),L=1,6),FR(N,1), * FR(N,3),ERROR(n),'--excl' endif else c c...deposition (optional) if(iexdep.eq.1)then ndig=10 if(error(n).gt.0.01d0)then write(cdum1,'(f10.3)')fr(n,3) write(cdum2,'(f10.3)')error(n) else write(cdum1,'(f10.4)')fr(n,3) write(cdum2,'(f10.4)')error(n) endif call lzero(cdum1,ndig) call lzero(cdum2,ndig) if(fr(n,4).eq.0.d0)then if(error(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'---excl ' else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'--excl ' endif else if(error(n).gt.0.01d0)then WRITE(NOUT,68)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'---excl ',fr(n,4) else WRITE(NOUT,168)(JK(N,L),L=1,6),FR(N,1),cdum1(1:10), * cdum2(1:10),'--excl ',fr(n,4) endif endif endif c endif endif lcount=lcount+1 if(NOUT.eq.0.and.lcount.eq.lmax)goto 700 goto 7 c 700 write(*,716) read(*,'(i1)',err=705)jj if(jj.eq.1)goto 709 705 lcount=0 c 7 CONTINUE if(NOUT.eq.0.and.lcount.gt.lmax-12)then write(*,'('' Press ENTER to continue '',$)') read(*,'(i1)',err=709)jj endif c 709 if(k.ne.5)then if(nout.eq.0)then WRITE(NOUT,53) else WRITE(NOUT,553) endif endif c 81 FORMAT(I4,'.',A,2(I5,2I4,1x),F15.4,F14.4,F9.4,2a) 181 FORMAT(I4,'.',A,2(I5,2I4,1x),F15.4,' B',F12.4,F9.4,a, * 1PE8.1,0PF14.4) 721 FORMAT(I5,'.',A,2(I5,2I4,1x),F15.4,F14.4,F9.4,2a) 68 FORMAT(1X,2(i3,','),i3,' <-',2(I3,','),i3,F17.3,3a,1PE10.2) 168 FORMAT(1X,2(i3,','),i3,' <-',2(I3,','),i3,F18.4,3a,1PE9.2) IF(NOUT.NE.0.and.k.ne.5)WRITE(NOUT,410) * SDEV,NINFIT,SWDEV,NINFIT-K1 410 FORMAT(/' Sigma = ',F14.6,28X,I4,' Distinct lines in fit'/ 1 ' Sigma.w = ',F14.6,28X,I4,' Degrees of freedom') C C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...CONSTANTS C C...On output to disk additional conversion of constants into customary units C and writing in the convention of error given in units of last quoted C digit in constant value (three digit error) C IF(NOUT.EQ.3.and.k.ne.5)THEN WRITE(NOUT,'(1X)') L=0 DO 1002 N=1,NCONST UNIT=NAMUN(UNITS(N)) ASCLD=A(N)*10.D0**(3*(UNITS(N)-1)) IF(IA(N).EQ.0)GOTO 1003 C C...value with error L=L+1 ESCLD=C(L)*10.D0**(3*(UNITS(N)-1)) WRITE(CONVAL,'(F25.16)')ASCLD WRITE(ERVAL,'(F25.16)')ESCLD CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR) if(ndcon+nderor.gt.22)ndcon=22-nderor WRITE(OUTCON(N),411)T(N)(1:5),UNIT,CONVAL(1:NDCON), * ERVAL(1:NDEROR) GOTO 1002 C C...value without error, formatting is carried out to insert zero before C the decimal point and to discard digits from rounding errors 1003 IF(ASCLD.NE.0.D0)ASCLD=ASCLD+MOD(ASCLD,1.D-15*ASCLD) WRITE(CONVAL,'(F25.16)')ASCLD C C...correct leading zero IF(ASCLD.EQ.0.0D0)CONVAL=' 0.0 ' IF(CONVAL(8:8).EQ.' ')CONVAL(8:8)='0' IF(CONVAL(8:8).EQ.'-')THEN CONVAL(8:8)='0' CONVAL(7:7)='-' ENDIF C C...truncate to significant digits (locate a string of five zeroes) idot=9 iend=25 if(abs(ascld).lt.1.d0)then do 1105 m=10,25 if(conval(m:m).ne.'0')then mstart=m goto 1114 endif 1105 continue else mstart=idot+1 goto 1114 endif goto 1115 c 1114 do 1104 m=mstart,21 if(conval(m:m+4).eq.'00000')then iend=m-1 goto 1115 endif 1104 continue c 1115 WRITE(OUTCON(N),1005)T(N)(1:5),UNIT,CONVAL(1:iend) 1002 CONTINUE c LLLL=NCONST/2 DO 1107 N=1,LLLL WRITE(NOUT,1108)OUTCON(N),OUTCON(N+LLLL) 1107 CONTINUE IF(2*LLLL.NE.NCONST)WRITE(NOUT,1108)' ',OUTCON(NCONST) ENDIF 1108 FORMAT(1X,A38,1X,A38) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C IF(NOUT.EQ.0)GOTO 24 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C Variance covariance matrix (not printed) and correlation coefficients C (lower triangle only). Also mean values of correlation coefficients. C IF(K.LT.1.or.k.eq.5)GOTO 31 SDEVSQ=SWDEV*SWDEV DO 200 N=1,K1 DO 200 L=N,K1 H(N,L)=DD(N,L)*SDEVSQ 200 CONTINUE DO 201 N=1,K1 DO 201 L=1,N DENOM=C(L)*C(N) IF(DENOM.EQ.0.D0)DENOM=1.D-50 H(N,L)=H(L,N)/DENOM 201 CONTINUE C WRITE(NOUT,202) 202 FORMAT(/' CORRELATION COEFFICIENTS, C.ij:') DO 3035 N=1,20 NCC(N)=0 3035 CONTINUE NCORCO=0 CORCAV=0.D0 CCSUM=0.0D0 M=1 305 MM=M+7 IF(MM.GT.K1)MM=K1 WRITE(NOUT,301) WRITE(NOUT,230)(T(TPT(L)),L=M,MM) 230 FORMAT(7X,8(3X,A5,1X)) WRITE(NOUT,301) DO 302 N=M,K1 LLLL=MM IF(LLLL.GT.N)LLLL=N WRITE(NOUT,303)T(TPT(N)),(H(N,L),L=M,LLLL) DO 3030 L=M,LLLL IF(L.NE.N)THEN NCORCO=NCORCO+1 CORCAV=CORCAV+DABS(H(N,L)) CCSUM=CCSUM+H(N,L) NNN=DABS(H(N,L))/0.05+1 IF(NNN.EQ.21)NNN=20 NCC(NNN)=NCC(NNN)+1 IF(NCC(NNN).GT.67)NCC(NNN)=67 ENDIF 3030 CONTINUE 302 CONTINUE 303 FORMAT(1X,A5,1X,8(F9.4)) IF(MM.GE.K1)THEN GOTO 3031 ELSE M=M+8 GOTO 305 ENDIF 3031 IF(NCORCO.NE.0)THEN WRITE(NOUT,3032)CORCAV/NCORCO,CCSUM/NCORCO 3032 FORMAT(1X/' Mean value of |C.ij|, i.ne.j =', * F9.4/ * ' Mean value of C.ij, i.ne.j =', * F9.4/) ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C...bar graph of distribution of correlation coefficients (disk only) IF(K.eq.2.or.k.eq.3)THEN WRITE(NOUT,'(1X/'' Distribution of values of |C.ij|, i.ne.j:'' * //10X,6(8X,I2)/10X,6(''....,....|''),''....,..'' * )')(N,N=10,60,10) DO 3036 N=1,20 WRITE(NOUT,3037)NCC(N),N*0.05,('*',LLLL=1,NCC(N)) 3036 CONTINUE 3037 FORMAT(1X,I2,' <',F4.2,1X,67A1) ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...MOST INFLUENTIAL TRANSITIONS and WORST FITTING LINES C C Most infl. trans. are defined by magnitude of diagonal elements hjj of C H = J (J'J)**(-1) J' C C here the design (derivative) matrix J is D, and (J'J)**(-1) is DD, C in ZK thesis J is X C IF(NOUT.EQ.3.and.k.ne.5)THEN call hjjset(K1,hjjav) c hjjcut=0.5 write(nout,'(79(1h_))') IF(work(i).lt.hjjcut)then write(nout,609)hjjcut goto 608 else write(nout,610)hjjcut endif 609 format(1x/' There are no influential transitions with hjj >', * f6.3) 610 FORMAT(1x/' MOST INFLUENTIAL TRANSITIONS:'/1x, * T42,'Observed',T56,'Obs-Calc',T70,'hjj >',F5.3/) do 607 L=I,1,-1 LL=ipt(L) if(linfit(LL).eq.0)goto 607 if(work(L).lt.hjjcut)goto 608 write(nout,606)LL,(JK(LL,LLL),LLL=1,6),FR(LL,1),FR(LL,3), * hjj(ll) 607 continue 606 FORMAT(I4,'. ',2(I5,2I4,1x),F15.4,F14.4,F9.4) 608 write(nout,611)hjjav 611 format(1x/1x,T46,'mean of all values =',f7.4) c c...worst fitting lines write(NOUT,105) do 106 L=1,I if(LINFIT(L).gt.0)then work(L)=abs(fr(L,3)/error(L)) if(L.ne.1)then if( (linfit(L).eq.2.and.linfit(L-1).eq.2) .and. * (fr(L,1).eq.fr(L-1,1)) )work(L)=0.d0 endif else WORK(L)=0.d0 endif ipt(L)=L 106 continue int21=1 call sortc(int21,I) nhj=I-6 if(nhj.lt.1)nhj=1 write(NOUT,110)(ipt(L),FR(ipt(L),3)/error(ipt(L)),L=I,nhj,-1) 105 FORMAT(1x/' WORST FITTING LINES (line number=(o-c)/error):') 110 FORMAT(1x,7(I5,'=',F5.1)) c ENDIF C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C Contributions to frequencies (disk only) or hjj/DFBETAS values for C transitions C 304 IF(K-2.LT.1)GOTO 31 C C...heading and options C 501 write(*,500) 500 format(1x/ * 8x,'-1 = values of hjj/DFBETAS for transitions'/ * 8x,' 0 = All contributions to transition frequencies'/ * 8x,' 1 = Only contributions from selected constants'// * 8x,'.... ',$) read(*,'(i5)',err=501)iconse c c...hjj/DFBETAS - for exposition see: C C J.Demaison, J.Cosleou, R.Bocquet, and A.G.Lesarri, C J.Mol.Spectrosc. 167,400-418(1994) c if(iconse.eq.-1)then c c..Transitions and HJJ's c call hjjset(K1,hjjav) write(nout,'(79(1h_))') WRITE(NOUT,350) nspace=1 ncomm=2 write(*,'(1x/7x,'' '',$)') do 351 n=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 355 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 355 endif C if(linfit(N).eq.0)goto 351 write(NOUT,352)n,(jk(n,l),l=1,6),fr(n,1),fr(n,3), * error(n),fr(n,3)/error(n),hjj(n) 351 continue 350 FORMAT(/' FITTED TRANSITIONS, SCALED DEVIATIONS, HJJ''S:'// * 1x,t34,'obs',t45,'obs-calc',t56,'error ', * '(o-c)/Err',t76,'hjj'/) 352 FORMAT(i4,'.',2x,3i3,i4,2I3,F14.4,F12.4,f8.4,F9.3,F10.5) c C...HJJ's and DFBETAS C C The mode of DFBETAS calculation used here i.e. repetition of the whole C least squares fit for each line in turn is probably the slowest possible, C but it is a rare operation and is fast enough anyway - the Demaison C paper cites references, which probably contain a much faster route. C nspace=1 ncomm=2 write(nout,370) write(nout,376)(T(TPT(L))(1:5),L=1,K1) write(nout,'(1x)') do 371 n=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 375 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 375 endif c if(linfit(N).eq.0)goto 371 WRITE(*,802) 802 FORMAT(1H+,'.',$) CALL DFBET(DFBETA,K1,N) KFIN=K1 IF(KFIN.GT.11)KFIN=11 write(NOUT,373)n,hjj(n),(dfbeta(L),L=1,KFIN) 377 IF(KFIN.LT.K1)THEN KST=KFIN+1 KFIN=KST+10 IF(KFIN.GT.K1)KFIN=K1 write(NOUT,374)(dfbeta(L),L=KST,KFIN) GOTO 377 ENDIF 371 continue write(*,'(1x)') goto 31 endif 370 FORMAT(/' HJJ''S AND DFBETAS:'/) 376 FORMAT(8x,'hjj ',11(1x,a)/13x,11(1x,a)/13x,11(1x,a)/ * 13x,11(1x,a)/) 373 FORMAT(i4,'.',F8.5,11F6.2) 374 FORMAT(13x,11F6.2) c c...continuation of selection for contributions c if(iconse.lt.0.or.iconse.gt.1)goto 501 if(iconse.eq.1)then 502 write(*,503) 503 format(1x/' Please specify up to 10 comma separated numbers' * ' of constants for contributions'/' (smaller number of constants' * ,' is to be padded out with commas)'/10x,' ',$) do 515 n=1,10 namcd(n)=0 515 continue read(*,*,err=502)namcd noutc=0 do 504 n=1,10 if(namcd(n).ne.0)then noutc=noutc+1 if(namcd(n).lt.0.or.namcd(n).gt.nconst. * or.ia(namcd(n)).ne.1)namcd(n)=tpt(1) else goto 505 endif 504 continue 505 write(*,506)(t(namcd(n))(1:5),n=1,noutc) 506 format(1x/' Derivatives will be output for following', * ' constants:'/11x,10a) 508 write(*,507) 507 format(30x,' OK (1/0) ? ',$) read(*,'(i5)',err=508)iok if(iok.ne.1)goto 502 endif C C...CONTRIBUTIONS - the F descriptor in output format C is chosen individually for each derivative according to its magnitude C C...values of contributions DO 203 N=1,I DO 203 L=1,K1 D(N,L)=D(N,L)*A(TPT(L)) 203 CONTINUE c write(nout,'(79(1h_))') WRITE(NOUT,204) 204 FORMAT(/' CONTRIBUTIONS OF INDIVIDUAL CONSTANTS TO FREQUENCY:'/) if(iconse.eq.1)then write(nout,509)(t(namcd(n))(1:5),n=1,noutc) 509 format(1x,t32,'obs',t45,'o-c IFIT',2x,10(3x,a,3x)) write(nout,'(1x)') do 511 n=1,i write(NOUT,510)n,(jk(n,l),l=1,6),fr(n,1),fr(n,3),linfit(n), * (d(n,tptinv(namcd(l))),l=1,noutc) 511 continue goto 514 endif 510 FORMAT(i4,'.',2x,3i3,i4,2I3,2F12.3,i2,' ',10F11.3) c IST=1 ISP=8 390 IF(ISP.GT.K1)ISP=K1 WRITE(NOUT,230)(T(TPT(L)),L=IST,ISP) IF(ISP.LT.K1)THEN IST=IST+8 ISP=ISP+8 GOTO 390 ENDIF WRITE(NOUT,301) 301 FORMAT(1X) nspace=1 ncomm=2 DO 320 N=1,I if(ispace(nspace).eq.n)then nspace=nspace+1 write(NOUT,'(1x)') endif c 380 if(icomml(ncomm).eq.n)then write(NOUT,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 380 endif C FORMAL=F1 IST=1 315 ISP=IST+7 IF(ISP.GT.K1)ISP=K1 314 DO 310 L=IST,ISP C C...determine LLL in format descriptor F9.LLL such that largest number C of significant digits are displayed for the current contribution LLL=5 DER=ABS(D(N,L)) DIV=10. 312 IF(DER/DIV.LT.1.)GOTO 311 DIV=DIV*10. LLL=LLL-1 GOTO 312 311 IF(LLL.LT.0)LLL=0 LLLL=L IF(LLLL.GT.8)LLLL=L-IST+1 FORMAR(13+LLLL*5)=CHAR(LLL+ICHAR('0')) 310 CONTINUE IF(ISP.LE.8)THEN WRITE(NOUT,FORMAL)N,(D(N,L),L=IST,ISP) ELSE WRITE(NOUT,FORMAL)(D(N,L),L=IST,ISP) ENDIF IF(ISP.EQ.K1)GOTO 320 FORMAL=F2 IST=IST+8 GOTO 315 C 320 CONTINUE C C...restore original derivatives C 514 DO 512 N=1,I DO 512 L=1,K1 D(N,L)=D(N,L)/A(TPT(L)) 512 CONTINUE c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c c...end of results output c 31 if(k.eq.5)then write(nout,'(79(1H_)/79(1H_))') else if(k.le.3)write(nout,'(79(1h_))') endif NOUT=0 CLOSE(3) GOTO 23 c 30 WRITE(NOUT,405)SDEV,SWDEV 24 IF(IOUT.NE.1.AND.IOUT.NE.2)GOTO 300 IOUT=0 GOTO 16 C c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c C O U T P U T M E N U C 300 WRITE(*,22) 22 FORMAT(/' OUTPUT OF RESULTS - CHOOSE ONE OF:'// 1 9X,'0 = no output'/ 1 9X,'1 = standard output'/ 1 9X,'2 = above + distribution plot of corr. coeff.'/ 1 9X,'3 = above + freq. contributions or hjj/DFBETAS'/ 1 9X,'4 = constants file for ASROT'/ 1 9X,'5 = frequency listing in deposition format'/ 1 9x,'6 = generate .LIN and .PAR files for SPFIT .... ',$) READ(*,15,ERR=300)K IF(K.LT.1.OR.K.GT.6)GOTO 23 NOUT=3 1901 WRITE(*,17) if(K.ne.6)then READ(*,'(A)',ERR=1901)FILNAM OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1901) endif if(k.eq.5)then 401 write(*,400) 400 FORMAT(/' ARE EXCLUDED LINES ALSO TO BE LISTED (1/0): ',$) read(*,'(I5)',err=401)iexdep if(iexdep.lt.0.or.iexdep.gt.1)goto 401 endif C c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c c...output of a constants file C IF(K.EQ.4)THEN DO 6869 LL=NCONST,1,-1 IF(A(LL).NE.0)GOTO 6870 6869 CONTINUE LL=NCONST 6870 WRITE(3,6866)LL LLL=0 DO 6867 L=1,LL ERCON=0.D0 IF(IA(L).EQ.1)THEN LLL=LLL+1 ERCON=C(LLL) ENDIF WRITE(3,6868)T(L),A(L),ERCON 6867 CONTINUE 6866 FORMAT(' NCON =',I5) 6868 FORMAT(1X,A6,2F30.19) WRITE(3,'('' |''/3x,22(1H-),4A)')' Reduction-',NAMRED, * ', Representation-',NAMREP GOTO 31 ENDIF C c_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ c c...output of a .LIN file for SPFIT C IF(K.EQ.6)THEN 1401 write(*,1400) 1400 FORMAT(/' Generic name for .LIN and .PAR files: ',$) READ(*,'(A)',ERR=1401)FILNAM filnam=filnam(1:lentrm(filnam))//'.lin' OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1401) nspace=1 ncomm=2 ivnum=-1 7005 write(*,7004) 7004 format(1x/ ' Add a vibrational quantum number (0/1) ? ',$) vv='00' read(*,'(i5)',err=7005)kk if(kk.eq.1)then 7007 write(*,7006) 7006 format( ' Value of the vibrational quantum number = ',$) read(*,'(i5)')m if(m.lt.0.or.m.gt.9)goto 7007 ivnum=m write(vv,'(2i1)')ivnum,ivnum endif do 7000 kk=1,I linout='' if(ispace(nspace).eq.kk)then nspace=nspace+1 linout(1:2)=' !' endif c 7001 if(icomml(ncomm).eq.kk)then write(lincom,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 if(lentrm(linout).gt.0)then linout=linout(1:lentrm(linout))//lincom(1:lentrm(lincom)) else linout=' '//lincom(1:lentrm(lincom)) endif goto 7001 endif c if(fr(kk,4).ne.0.d0)then c write(linrel,'(f10.5)')fr(kk,4) write(linrel,'(1PE10.2)')fr(kk,4) else linrel=' 1.00000' endif lendes=lentrm(linout) if(lendes.eq.1.and.linout(l:1).le.' ')lendes=0 lendsc=lentrm(descr(kk)) if(lendsc.eq.1.and.descr(kk)(l:1).le.' ')lendsc=0 if(lendsc.gt.0)then if(lendes.eq.0)then linout=' #'//descr(kk)(1:lendsc) else linout=linout(1:lendes)//' #'//descr(kk)(1:lendsc) endif lendes=lentrm(linout) endif c errprn=error(kk) if(linfit(kk).eq.0)errprn=error(kk)+900.d0 if(lendes.eq.0)then if(ivnum.lt.0)WRITE(3,7002)(JK(kK,M),M=1,6),FR(kK,1),ERRprn, * linrel if(ivnum.ge.0)WRITE(3,7003)(JK(kK,M),M=1,3),ivnum, * (JK(kk,m),m=4,6),ivnum,FR(kK,1),ERRprn, * linrel else if(ivnum.lt.0)WRITE(3,7002)(JK(kK,M),M=1,6),FR(kK,1),ERRprn, * linrel,linout(1:lendes) if(ivnum.ge.0)WRITE(3,7003)(JK(kK,M),M=1,3),ivnum, * (jk(kk,m),m=4,6),ivnum,FR(kK,1),ERRprn, * linrel,linout(1:lendes) endif 7000 continue 7002 format(6i3,' 0 0 0 0 0 0',f15.5,f12.5,2a) 7003 format(8i3,' 0 0 0 0',f15.5,f12.5,2a) close(3) c c...output of a .PAR file for SPFIT C k=lentrm(filnam) filnam=filnam(1:k-4)//'.par' OPEN(3,FILE=FILNAM,STATUS='UNKNOWN',ERR=1901) c write(3,'(72a1)')com do 1405 ncon=nconst,1,-1 if(a(ncon).ne.0.d0)goto 1406 1405 continue 1406 write(3,1407)ncon,i,5,0, * ' 0.0000E+000 1.0000E+002 1.0000E+000 1.0000000000' 1407 format(i5,i6,2i5,a) if(namred.eq.'A,i2,a')then if(namrep.eq.'I.r ')then write(3,'(a,i2,a)') * 'a 1 ',ivnum+1,' 0, , , , , , , ,,,,,' else write(3,'(a,i2,a)') * 'a 1 ',-(ivnum+1),' 0, , , , , , , ,,,,,' endif else if(namrep.eq.'I.r ')then write(3,'(a,i2,a)') * 's 1 ',ivnum+1,' 0, , , , , , , ,,,,,' else write(3,'(a,i2,a)') * 's 1 ',-(ivnum+1),' 0, , , , , , , ,,,,,' endif endif c do 1408 kk=1,ncon picon(kk)=a(kk) if(kk.ge.4.and.kk.le.8)picon(kk)=-a(kk) if((kk.eq.7.or.kk.eq.8).and.namred.eq.'S')picon(kk)=a(kk) if(namrep.eq.'III.r')then if(kk.eq. 7)picon(kk)=-picon(kk) if(kk.eq. 8)picon(kk)=-picon(kk) if(kk.eq.13)picon(kk)=-picon(kk) if(kk.eq.14)picon(kk)=-picon(kk) if(kk.eq.15)picon(kk)=-picon(kk) if(kk.eq.21)picon(kk)=-picon(kk) if(kk.eq.22)picon(kk)=-picon(kk) if(kk.eq.23)picon(kk)=-picon(kk) if(kk.eq.24)picon(kk)=-picon(kk) if(kk.eq.31)picon(kk)=-picon(kk) if(kk.eq.32)picon(kk)=-picon(kk) if(kk.eq.33)picon(kk)=-picon(kk) if(kk.eq.34)picon(kk)=-picon(kk) if(kk.eq.35)picon(kk)=-picon(kk) endif if(namred.eq.'S'.and.kk.eq. 8)pinums(kk)='500' if(namred.eq.'S'.and.kk.eq.14)pinums(kk)='501' if(namred.eq.'S'.and.kk.eq.15)pinums(kk)='600' if(namred.eq.'S'.and.kk.eq.22)pinums(kk)='502' if(namred.eq.'S'.and.kk.eq.23)pinums(kk)='601' if(namred.eq.'S'.and.kk.eq.24)pinums(kk)='700' if(namred.eq.'S'.and.kk.eq.32)pinums(kk)='503' if(namred.eq.'S'.and.kk.eq.33)pinums(kk)='602' if(namred.eq.'S'.and.kk.eq.34)pinums(kk)='701' if(namred.eq.'S'.and.kk.eq.35)pinums(kk)='800' 1408 continue c do 1402 kk=1,ncon if(a(kk).eq.0.d0.and.ia(kk).eq.0)goto 1402 itstrt=1 if(t(kk)(1:1).eq.' ')itstrt=2 if(ia(kk).eq.0)then parwt=1.D-35 else parwt=1.D+20 endif if(kk.ge.4.and.kk.le.8)then if((kk.eq.7.or.kk.eq.8).and.namred.eq.'S')then write(3,1403)pinums(kk)//vv,picon(kk),parwt,t(kk)(itstrt:5) else write(3,1404)pinums(kk)//vv,picon(kk),parwt,t(kk)(itstrt:5) endif else write(3,1403)pinums(kk)//vv,picon(kk),parwt,t(kk)(itstrt:5) endif 1402 continue write(3,1410) close(3) NOUT=0 c 1403 format(11x,a5,1PE24.15,E16.8,' /',a) 1404 format(11x,a5,1PE24.15,E16.8,' /-',a) 1410 format('!'/'!',78(1H-)/'!'/ * '! Generated automatically by ASFIT'/'!'/ * '!',78(1H-)/) c write(*,'(1x/'' SUCCESSFUL OUTPUT OF: '',a/18x,''AND: '',a/)') * filnam(1:k-4)//'.lin',filnam(1:lentrm(filnam)) GOTO 23 ENDIF c GOTO 25 C C_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ C C...Output of an .ASF data file C 23 WRITE(*,14) 14 FORMAT(/' DATA FILE TO BE WRITTEN TO DISK (1/0) ? ',$) READ(*,15,ERR=23)L 15 FORMAT(I5) IOUT=1 IF(L.NE.1)GOTO 24 100 WRITE(*,17) 17 FORMAT(/' NAME TO BE USED FOR THE OUTPUT FILE : ',$) READ(*,'(A)',ERR=100)FILNAM OPEN(2,FILE=FILNAM,STATUS='UNKNOWN',ERR=100) WRITE(2,51)COM write(2,'(22(1H-),5(1H+),61(1H-))') do 101 ncon=nconst,1,-1 if(a(ncon).ne.0.d0)goto 102 101 continue 102 WRITE(2,54)'number of transitions:',I write(2,54)'number of constants :',ncon write(2,'(5(1H-),30(1H+),53(1H-))') do 103 k=1,ncon if(k.eq.1)then WRITE(2,150)IA(K),A(K),T(K)(1:5),' Reduction-', * namred,', Representation-',namrep else WRITE(2,150)IA(K),A(K),T(K)(1:5) endif 103 continue c write(2,55) write(2,56) write(2,554) nspace=1 ncomm=2 do 3400 kk=1,I if(ispace(nspace).eq.kk)then nspace=nspace+1 write(2,'(1x)') endif c 133 if(icomml(ncomm).eq.kk)then write(2,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 133 endif c lendes=lentrm(descr(kk)) if(lendes.eq.0)then if(fr(kk,4).eq.0.d0)then WRITE(2,140)(JK(kK,M),M=1,6),FR(kK,1),ERROR(kK),LINFIT(kK) else WRITE(2,140)(JK(kK,M),M=1,6),FR(kK,1),ERROR(kK),LINFIT(kK), * FR(kk,4) endif else if(fr(kk,4).eq.0.d0)then WRITE(2,141)(JK(kK,M),M=1,6),FR(kK,1),ERROR(kK),LINFIT(kK), * DESCR(kK)(1:lendes) else WRITE(2,140)(JK(kK,M),M=1,6),FR(kK,1),ERROR(kK),LINFIT(kK), * fr(kk,4),DESCR(kK)(1:lendes) endif endif 3400 continue write(2,554) write(2,55) CLOSE(2) c 51 FORMAT(1X,72A1) 53 format(79(1H-)) 553 format(102(1H-)) 554 format(3(5(1H+),5(1H-)),20(1H+),15(1H-),5(1H+),10(1H-),8(1H+)) 54 format(a22,i5) 55 format(88(1H-)) 56 format(3x,'J'' K-1'' K+1'' J K-1 K+1 frequency', * ' error fit? blend Note') 140 FORMAT(6I5,F20.6,F15.6,I5,1PE9.2,1x,A) c140 FORMAT(6I5,F20.6,F15.6,I5,F9.4,1x,A) 141 FORMAT(6I5,F20.6,F15.6,I5,10x,A) 150 FORMAT(I5,F30.19,3X,5A) C IOUT=1 GOTO 24 C 16 RETURN END C C_____________________________________________________________________________ C SUBROUTINE DFBET(Z,K1,N) C C DFBETAS for K1 fitted constants associated with the N'th transition C (the K1 calculated values are passed back in array Z) C C Changes in constants relative to last iteration are in B, their errors C are in C, D'Y is in K1+1'st column of H, (D'D) is set up in top-left C corner of H and inv(D'D) placed in DD C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250, MAXDIM=JMAX/2+1, * NNJA=2*NLINES+1) COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED COMMON /BLKH/H * /BLKR/R * /BLK2/NJA,B * /BLK3/C,D,DD,HJJ INTEGER*2 JK(NLINES,6),IA(NCONST),LINFIT(NLINES), * NJA(NNJA,4) REAL*8 FR(NLINES,4),A(NCONST),C(NCONST),H(MAXDIM,MAXDIM), * R(MAXDIM,MAXDIM),DD(NCONST,NCONST), * Z(NCONST),B(NCONST),BNEW(NCONST) REAL*4 D(NLINES,NCONST),ERROR(NLINES),WEIGHT(NLINES), * HJJ(NLINES) character*8 descr(nlines) C if(linfit(N).eq.0)goto 300 C C...Evaluate inv(D'D) from the old derivative matrix in D, but without C line N - inversion is through diagonalisation as in M/P C C array Z is used: C C 1. to hold scaling constants improving conditioning in evaluation of C inv(D'D) C 2. for DFBETAS associated with line N C DO 30 J=1,K1 DO 28 K=J,K1 E=0.D0 DO 27 L=1,I IF(L.eq.N)GOTO 27 IF(LINFIT(L).eq.0)GOTO 27 E=E+D(L,J)*D(L,K)*WEIGHT(L) 27 CONTINUE H(J,K)=E 28 CONTINUE Z(J)=SQRT(1.D0/H(J,J)) E=0.D0 DO 29 L=1,I IF(L.eq.N)GOTO 29 IF(LINFIT(L).eq.0)GOTO 29 E=E+D(L,J)*FR(L,3)*WEIGHT(L) 29 CONTINUE H(J,K1+1)=E 30 CONTINUE C DO 301 J=1,K1 DO 301 K=J,K1 H(J,K)=H(J,K)*Z(J)*Z(K) H(K,J)=H(J,K) 301 CONTINUE IEGEN=0 C CALL HDIAG(K1,IEGEN,L3) C DO 35 J=1,K1 BNEW(J)=0.D0 DO 34 K=1,K1 E=0.D0 DO 33 L=1,K1 E=E+R(J,L)*R(K,L)/H(L,L) 33 CONTINUE E=E*Z(J)*Z(K) DD(J,K)=E BNEW(J)=BNEW(J)+E*H(K,K1+1) 34 CONTINUE 35 CONTINUE C C...DFBETAS (from difference of corrections B and BNEW to constants from C previous least squares iteration) C 300 if(LINFIT(N).eq.0)then do 360 L=1,K1 Z(L)=0.0d0 360 continue return else do 361 L=1,K1 Z(L)=( B(L)-BNEW(L) )/C(L) 361 CONTINUE endif C RETURN END C_____________________________________________________________________________ C subroutine hjjset(K1,hjjav) C C Calculation of the HJJ coefficients describing the influence of transition C J on the fit: C transitions with hjj > 0.5 should be used with care C transitions with hjj > 0.9 strongly skew the fit and should really only C be tolerated in intermediate stages of constant determination C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NLINES=10000, NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) COMMON /BLK1/JK,FR,ERROR,WEIGHT,LINFIT,A,IA,DESCR,I,IOUT,IRED * /BLK3/C,D,DD,HJJ * /SORTCC/WORK,IPT INTEGER*2 JK(NLINES,6),IA(NCONST),LINFIT(NLINES),IPT(NLINES) REAL*4 D(NLINES,NCONST),ERROR(NLINES),WEIGHT(NLINES), * HJJ(NLINES) REAL*8 FR(NLINES,4),A(NCONST),C(NCONST), * DD(NCONST,NCONST),WORK(NLINES) character*8 descr(nlines) C C...set up hjj's = diagonal elements of H = D inv(D'D) D' C do 600 L=1,I if(linfit(L).eq.0)goto 600 do 601 LL=1,K1 work(LL)=0.0 do 602 LLL=1,K1 work(LL)=work(LL)+D(L,LLL)*DD(LLL,LL) 602 continue 601 continue hjj(L)=0.0 do 603 LL=1,K1 hjj(L)=hjj(L)+work(LL)*D(L,LL) 603 continue hjj(L)=hjj(L)*weight(L) 600 continue C C...sort hjj's according to magnitude C hjjcut=0.0 hjjav=0.0 nhjj=0 do 605 L=1,I if(linfit(L).gt.0)then work(L)=hjj(L) hjjav=hjjav+hjj(L) nhjj=nhjj+1 else work(L)=0.0d0 endif ipt(L)=L 605 continue if(nhjj.gt.0)hjjav=hjjav/nhjj int21=1 call sortc(int21,i) c return end C C_____________________________________________________________________________ C subroutine lzero(number,ndig) character*20 number integer*2 n,ndig c c Addition of a leading zero to a number in a string c do 1 n=1,ndig-2 if(number(n:n+1).eq.' .')number(n:n+1)='0.' if(number(n:n+2).eq.' -.')number(n:n+2)='-0.' 1 continue c return end C_____________________________________________________________________________ C SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR) C C Constant and error formating for output C IMPLICIT INTEGER*2 (I-N) CHARACTER*25 CONVAL,ERVAL,CONOUT,EROUT C NERD=2 NDCON=0 NDEROR=0 NDNOTZ=0 ICZERO=0 C C...Go through digits of constant value adding those to output buffer while C checking at the same time digits of the error value, reacting as C necessary. Terminate when either three digits in the error value are C transferred or, if error has more than three digits before the decimal C point, the decimal point is reached. DO 1 N=1,25 NDCON=NDCON+1 CONOUT(NDCON:NDCON)=CONVAL(N:N) C C...ensure that zero precedes the decimal point and use ICZERO to monitor C whether decimal point has been reached IF(CONVAL(N:N).EQ.'.')ICZERO=1 if(ndcon.gt.1.and.n.gt.1)then IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.' ') * CONOUT(NDCON-1:NDCON-1)='0' endif if(ndcon.gt.2.and.n.gt.1)then IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.'-')THEN CONOUT(NDCON-2:NDCON-2)='-' CONOUT(NDCON-1:NDCON-1)='0' ENDIF endif C C...use NDNOTZ (number of digits not zero) to monitor whether significant C digits in constant value have been reached IF(CONVAL(N:N).GE.'1'.AND.CONVAL(N:N).LE.'9') * NDNOTZ=NDNOTZ+1 C C...do not transfer error digit if it is a leading zero, dot or space IF(NDEROR.EQ.0 .AND. (ERVAL(N:N).EQ.' '.OR. * ERVAL(N:N).EQ.'0'.OR.ERVAL(N:N).EQ.'.') )GOTO 1 C C...exit if error larger than value and decimal point reached IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 C C...do not transfer the dot in error value IF(ERVAL(N:N).EQ.'.')GOTO 1 C C...transfer error digits until three or, if more, the first significant C digit in value is reached NDEROR=NDEROR+1 EROUT(NDEROR:NDEROR)=ERVAL(N:N) IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 1 CONTINUE 2 CONVAL(1:NDCON)=CONOUT(1:NDCON) if(nderor.lt.6.and.nderor.ge.1)then ERVAL(1:NDEROR)=EROUT(1:NDEROR) else if(nderor.eq.0)then nderor=1 erval(1:1)='*' endif if(nderor.ge.6)then nderor=5 ERVAL(1:5)='*****' endif endif C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ASREDA(I,J,N,A,IRED) C C Evaluation of a Wang submatrix for reduction A, where: C C I - on output contains the dimensionality of the calculated subm. C J - J quantum number for the subm. C N - number of the required subm. = 1,2,3,4 for E+,E-,O+,O- resp. C A - array containing rotational and cd. constants in usual order C C Evaluated submatrix is placed in H C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) INTEGER*4 JJ1,K,KK COMMON /BLKH/H * /BLKR/U DIMENSION A(NCONST),H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM) C JJ1=J JJ1=JJ1*(JJ1+1) Q=JJ1 QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C IF(IRED.GT.0)THEN C=0.5D0*(A(2)+A(3)) DK2=A(1)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ELSE C=0.5D0*(A(2)+A(1)) DK2=A(3)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ENDIF DK0=C*Q-A(4)*QQ+A(9)*QQQ+A(16)*QQQQ+A(25)*QQQQQ DK4=-A(6)+A(11)*Q+A(18)*QQ+A(27)*QQQ DK6=A(12)+A(19)*Q+A(28)*QQ DK8=A(20)+A(29)*Q DK10=A(30) C IF(IRED.GT.0)THEN OK0=0.25D0*(A(2)-A(3))-A(7)*Q+A(13)*QQ+A(21)*QQQ+ * A(31)*QQQQ ELSE OK0=0.25D0*(A(1)-A(2))-A(7)*Q+A(13)*QQ+A(21)*QQQ+ * A(31)*QQQQ ENDIF OK2=0.5D0*(-A(8)+A(14)*Q+A(22)*QQ+A(32)*QQQ) OK4=0.5D0*(A(15)+A(23)*Q+A(33)*QQ) OK6=0.5D0*(A(24)+A(34)*Q) OK8=0.5D0*A(35) C C KA and KE are the lowest and highest K in the submatrix: C KA = 0, 2, 1, 1 for E+, E-, O+, O- resp. C K=(N-1)/2 M=N/2*2 KA=(4-N)/2*M+K KE=(J+K)/2*2-K M=(KE-KA)/2+1 C DO 1 I=1,M DO 1 K=1,M H(I,K)=0.D0 1 CONTINUE H(1,1)=0.D0 H(1,2)=0.D0 C C For summing over K loop variable KK equal to K+1 is used to avoid C loop count from zero for E+ C Variable I is a pointer to elements in the Wang submatrix C I=0 KA1=KA+1 KE1=KE+1 DO 3 KK=KA1,KE1,2 K=KK-1 I=I+1 FK2=K**2 H(I,I)=((((DK10*FK2+DK8)*FK2+DK6)*FK2+DK4)*FK2+DK2)*FK2 * +DK0 IF(K.LT.KE)GOTO 2 GOTO 3 2 C=(K+2)**2 FK8=FK2**4+C**4 FK6=FK2**3+C**3 FK4=FK2*FK2+C*C FK2=C+FK2 FLTB=JJ1-K*KK FLTC=JJ1-KK*(K+2) F2=FLTB*FLTC H(I,I+1)=(OK0+FK2*OK2+FK4*OK4+FK6*OK6+FK8*OK8)*DSQRT(F2) H(I+1,I)=H(I,I+1) 3 CONTINUE C C Final touch up: 4 -> O+ and O-, 6 -> E-, 5 -> E+ C IF(2-N)4,6,5 C C FLTA = +1 for O+, -1 for O- C 4 FLTA=-2*N+7 H(1,1)=FLTA*(OK0+2.D0*(OK2+OK4+OK6+OK8))*Q+H(1,1) GOTO 6 C 5 H(1,2)=H(1,2)*1.414213562373095D0 H(2,1)=H(1,2) 6 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ASREDS(I,J,N,A,IRED) C C Evaluation of a Wang submatrix for reduction S, where: C C I - on output contains the dimensionality of the calculated subm. C J - J quantum number for the subm. C N - number of the required subm. = 1,2,3,4 for E+,E-,O+,O- resp. C A - array containing rotational and cd. constants in usual order C C Evaluated submatrix is placed in H C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (NCONST=35, JMAX=250, MAXDIM=JMAX/2+1) INTEGER*4 K COMMON /BLKH/H * /BLKR/U DIMENSION A(NCONST),H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM) C Q=REAL(J+1)*REAL(J) QQ=Q*Q QQQ=QQ*Q QQQQ=QQQ*Q QQQQQ=QQQQ*Q C IF(IRED.GT.0)THEN C=0.5D0*(A(2)+A(3)) DK2=A(1)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ELSE C=0.5D0*(A(2)+A(1)) DK2=A(3)-C-A(5)*Q+A(10)*QQ+A(17)*QQQ+A(26)*QQQQ ENDIF DK0=C*Q-A(4)*QQ+A(9)*QQQ+A(16)*QQQQ+A(25)*QQQQQ DK4=-A(6)+A(11)*Q+A(18)*QQ+A(27)*QQQ DK6=A(12)+A(19)*Q+A(28)*QQ DK8=A(20)+A(29)*Q DK10=A(30) C IF(IRED.GT.0)THEN OKK2=0.25D0*(A(2)-A(3))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ ELSE OKK2=0.25D0*(A(1)-A(2))+A(7)*Q+A(13)*QQ+A(21)*QQQ+A(31)*QQQQ ENDIF OKK4=A(8)+A(14)*Q+A(22)*QQ+A(32)*QQQ OKK6=A(15)+A(23)*Q+A(33)*QQ OKK8=A(24)+A(34)*Q OKK10=A(35) C C KA and KE are the lowest and highest K in the submatrix: C KA = 0, 2, 1, 1 for E+, E-, O+, O- resp. C K=(N-1)/2 M=N/2*2 KA=(4-N)/2*M+K KE=(J+K)/2*2-K M=(KE-KA)/2+1 C DO 1 I=1,M DO 1 K=1,M H(I,K)=0.D0 1 CONTINUE H(1,1)=0.D0 H(1,2)=0.D0 C C For summing over K loop variable KK equal to K+1 is used to avoid C loop count from zero for E+ C Variable I is a pointer to elements in the Wang submatrix C I=0 KA1=KA+1 KE1=KE+1 DO 3 KK=KA1,KE1,2 K=KK-1 I=I+1 FK2=REAL(K)**2 c c... H(I,I)=((((DK10*FK2+DK8)*FK2+DK6)*FK2+DK4)*FK2+DK2)*FK2 * +DK0 c c... IF(K.GE.KE)GOTO 3 H(I,I+1)=OKK2*DSQRT(FPLUS(J,K,2)) H(I+1,I)=H(I,I+1) c c... IF(K.GE.KE-2)GOTO 3 H(I,I+2)=OKK4*DSQRT(FPLUS(J,K,4)) H(I+2,I)=H(I,I+2) c c... IF(K.GE.KE-4)GOTO 3 H(I,I+3)=OKK6*DSQRT(FPLUS(J,K,6)) H(I+3,I)=H(I,I+3) c c... IF(K.GE.KE-6)GOTO 3 H(I,I+4)=OKK8*DSQRT(FPLUS(J,K,8)) H(I+4,I)=H(I,I+4) c c... IF(K.GE.KE-8)GOTO 3 H(I,I+5)=OKK10*DSQRT(FPLUS(J,K,10)) H(I+5,I)=H(I,I+5) c 3 CONTINUE C C Final touch up: 4 -> O+ and O-, 12 -> E-, 5 -> E+ C IF(2-N)4,12,5 c c...O+ and O-, FLTA = +1 for O+, -1 for O- 4 FLTA=-2*N+7 H(1,1)=FLTA*OKK2*DSQRT(FMINUS(J,1,2))+H(1,1) IF(KE.LT.3)GOTO 6 H(1,2)=FLTA*OKK4*DSQRT(FMINUS(J,3,4))+H(1,2) H(2,1)=H(1,2) IF(KE.LT.5)GOTO 6 H(2,2)=FLTA*OKK6*DSQRT(FMINUS(J,3,6))+H(2,2) H(1,3)=FLTA*OKK6*DSQRT(FMINUS(J,5,6))+H(1,3) H(3,1)=H(1,3) IF(KE.LT.7)GOTO 6 H(1,4)=FLTA*OKK8*DSQRT(FMINUS(J,7,8))+H(1,4) H(4,1)=H(1,4) H(2,3)=FLTA*OKK8*DSQRT(FMINUS(J,5,8))+H(2,3) H(3,2)=H(2,3) IF(KE.LT.9)GOTO 6 H(1,5)=FLTA*OKK10*DSQRT(FMINUS(J,9,10))+H(1,5) H(5,1)=H(1,5) H(2,4)=FLTA*OKK10*DSQRT(FMINUS(J,7,10))+H(2,4) H(4,2)=H(2,4) H(3,3)=FLTA*OKK10*DSQRT(FMINUS(J,5,10))+H(3,3) GOTO 6 c c...root(2) elements of E+ 5 IF(KE.LT.2)GOTO 6 K=4 IF(KE.LT.6)K=1+KE/2 DO 10 KK=2,K H(1,KK)=H(1,KK)*1.414213562373095D0 H(KK,1)=H(1,KK) 10 CONTINUE c c...E- and remaining elements of E+, FLTA = +1 for E+ and -1 for E- 12 IF(KE.LT.2)GOTO 6 FLTA=-2*N+3 K=2 IF(N.EQ.2)K=1 H(K,K)= FLTA*OKK4*DSQRT(FMINUS(J,2,4))+H(K,K) IF(KE.LT.4)GOTO 6 H(K,K+1)=FLTA*OKK6*DSQRT(FMINUS(J,4,6))+H(K,K+1) H(K+1,K)=H(K,K+1) IF(KE.LT.6)GOTO 6 H(K,K+2)=FLTA*OKK8*DSQRT(FMINUS(J,6,8))+H(K,K+2) H(K+2,K)=H(K,K+2) H(K+1,K+1)=FLTA*OKK8*DSQRT(FMINUS(J,4,8))+H(K+1,K+1) IF(KE.LT.8)GOTO 6 H(K,K+3)=FLTA*OKK10*DSQRT(FMINUS(J,8,10))+H(K,K+3) H(K+3,K)=H(K,K+3) H(K+1,K+2)=FLTA*OKK10*DSQRT(FMINUS(J,6,10))+H(K+1,K+2) H(K+2,K+1)=H(K+1,K+2) 6 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION FPLUS(J,K,KDELTA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) INTEGER*4 K,KDELTA C C Values of f for matrix elements : C C dK-1 dK dK C + ___ ___ ___ C f (J,K) = | | (J-K-l) | | (J+K+l) = | | [J(J+1)-(K+l)(K+l-1)] C dK l=0 l=1 l=1 C FPLUS=1.0D0 RJ=J RK=K C DO 1 L=0,KDELTA-1 FPLUS=FPLUS*(RJ-RK-L) 1 CONTINUE DO 2 L=1,KDELTA FPLUS=FPLUS*(RJ+RK+L) 2 CONTINUE C RETURN END C C_____________________________________________________________________________ C FUNCTION FMINUS(J,K,KDELTA) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) INTEGER*4 K,KDELTA C C Values of f for matrix elements : C C dK-1 dK dK C - ___ ___ ___ C f (J,K) = | | (J+K-l) | | (J-K+l) = | | [J(J+1)-(K-l)(K-l+1)] C dK l=0 l=1 l=1 C FMINUS=1.0D0 RJ=J RK=K C DO 1 L=0,KDELTA-1 FMINUS=FMINUS*(RJ+RK-L) 1 CONTINUE DO 2 L=1,KDELTA FMINUS=FMINUS*(RJ-RK+L) 2 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HDIAG(N,IEGEN,NR) C C EIGENVALUES AND EIGENVECTORS OF A SQUARE, SYMMETRIC ARRAY C C N - dimensionality of array to be diagonalized C IEGEN - if =0 then calculate eigenvectors, if not then do not C NR - on output contains the number of iterations that were carried out C C Matrix to be diagonalized is placed in H C Eigenvalues are output on the diagonal of H, in ascending order C Eigenvectors are output in U in order of eigenvalues C C IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (JMAX=250, MAXDIM=JMAX/2+1) COMMON /BLKH/H * /BLKR/U DIMENSION H(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),X(MAXDIM),IQ(MAXDIM) C IF(IEGEN.NE.0)GOTO 15 DO 14 I=1,N DO 14 J=1,N IF(I-J.NE.0)GOTO 12 U(I,J)=1.D0 GOTO 14 12 U(I,J)=0.D0 14 CONTINUE 15 NR=0 IF(N-1.LE.0)GOTO 1000 NMI1=N-1 DO 30 I=1,NMI1 X(I)=0.D0 IPL1=I+1 DO 30 J=IPL1,N IF(X(I)-ABS(H(I,J)).GT.0.D0)GOTO 30 X(I)=ABS(H(I,J)) IQ(I)=J 30 CONTINUE 40 DO 70 I=1,NMI1 IF(I-1.LE.0)GOTO 60 IF(XMAX-X(I))60,70,70 60 XMAX=X(I) IPIV=I JPIV=IQ(I) 70 CONTINUE EPS=1.0D-12 IF(XMAX-EPS)1000,1000,148 148 CONTINUE NR=NR+1 130 IF(H(IPIV,IPIV)-H(JPIV,JPIV))140,141,142 140 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1 -SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.D0*H(IPIV,JPIV)**2)) GOTO 150 141 TANG=1.D0 GOTO 150 142 TANG=2.D0*H(IPIV,JPIV)/(H(IPIV,IPIV)-H(JPIV,JPIV) 1 +SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.D0*H(IPIV,JPIV)**2)) 150 COSINE=1.D0/SQRT(1.D0+TANG**2) SINE=TANG*COSINE HII=H(IPIV,IPIV) H(IPIV,IPIV)=COSINE**2*(HII+TANG*(2.D0*H(IPIV,JPIV)+TANG* 1 H(JPIV,JPIV))) H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.D0*H(IPIV,JPIV)- 1 TANG*HII)) H(IPIV,JPIV)=0.D0 IF(H(IPIV,IPIV)-H(JPIV,JPIV))153,153,152 152 HTEMP=H(IPIV,IPIV) H(IPIV,IPIV)=H(JPIV,JPIV) H(JPIV,JPIV)=HTEMP HTEMP=SIGN(1.D0,-SINE)*COSINE COSINE=ABS(SINE) SINE=HTEMP 153 CONTINUE DO 350 I=1,NMI1 IF(I-IPIV)210,350,200 200 IF(I-JPIV)210,350,210 210 IF(IQ(I)-IPIV.EQ.0)GOTO 240 IF(IQ(I)-JPIV)350,240,350 240 K=IQ(I) 250 HTEMP=H(I,K) H(I,K)=0.D0 IPL1=I+1 X(I)=0.D0 DO 320 J=IPL1,N IF(X(I)-ABS(H(I,J)).GT.0.D0)GOTO 320 X(I)=ABS(H(I,J)) IQ(I)=J 320 CONTINUE H(I,K)=HTEMP 350 CONTINUE X(IPIV)=0.D0 X(JPIV)=0.D0 DO 530 I=1,N IF(I-IPIV)370,530,420 370 HTEMP=H(I,IPIV) H(I,IPIV)=COSINE*HTEMP+SINE*H(I,JPIV) IF(X(I)-ABS(H(I,IPIV)).GE.0.D0)GOTO 390 X(I)=ABS(H(I,IPIV)) IQ(I)=IPIV 390 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF(X(I)-ABS(H(I,JPIV)).GE.0.D0)GOTO 530 X(I)=ABS(H(I,JPIV)) IQ(I)=JPIV GOTO 530 420 IF(I-JPIV)430,530,480 430 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(I,JPIV) IF(X(IPIV)-ABS(H(IPIV,I)).GE.0.D0)GOTO 450 X(IPIV)=ABS(H(IPIV,I)) IQ(IPIV)=I 450 H(I,JPIV)=-SINE*HTEMP+COSINE*H(I,JPIV) IF(X(IPIV)-ABS(H(I,JPIV)).GE.0.D0)GOTO 530 X(I)=ABS(H(I,JPIV)) IQ(I)=JPIV GOTO 530 480 HTEMP=H(IPIV,I) H(IPIV,I)=COSINE*HTEMP+SINE*H(JPIV,I) IF(X(IPIV)-ABS(H(IPIV,I)).GE.0.D0)GOTO 500 X(IPIV)=ABS(H(IPIV,I)) IQ(IPIV)=I 500 H(JPIV,I)=-SINE*HTEMP+COSINE*H(JPIV,I) IF(X(JPIV)-ABS(H(JPIV,I)).GE.0.D0)GOTO 530 510 X(JPIV)=ABS(H(JPIV,I)) IQ(JPIV)=I 530 CONTINUE IF(IEGEN.NE.0)GOTO 40 DO 550 I=1,N HTEMP=U(I,IPIV) U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV) U(I,JPIV)=-SINE*HTEMP+COSINE*U(I,JPIV) 550 CONTINUE GOTO 40 1000 CONTINUE C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HSORT(M) C C SORTING OF EIGENVALUES AND EIGENVECTORS COMING OUT OF HDIAG C - introduced after it was found that for representation III some C HDIAG eigenvalues did not keep ascending order. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*2 (I-N) PARAMETER (JMAX=250,MAXDIM=JMAX/2+1,NLINES=10000) C COMMON /BLKH/A * /BLKR/U * /SORTCC/WK,IPT REAL*8 A(MAXDIM,MAXDIM),U(MAXDIM,MAXDIM),WK(NLINES),D(MAXDIM) INTEGER*2 IPT(NLINES) isort=1 C C...Check whether additional sorting necessary C DO 150 I=1,M-1 IF(A(I+1,I+1).LT.A(I,I))GOTO 151 150 CONTINUE RETURN C C...Transfer eigenvalues from the diagonal of matrix A to D, and eigenvectors C from U to A C 151 DO 115 I=1,M D(I)=A(I,I) DO 115 J=1,M A(I,J)=U(I,J) 115 CONTINUE C C...Straightforward magnitude sorting C DO 122 I=1,M IPT(I)=I WK(I)=D(I) 122 CONTINUE CALL SORTC(isort,M) C C...Return of sorted eigenvectors in U and eigenvalues on the diagonal of A C 140 DO 117 II=1,M I=IPT(II) DO 117 J=1,M U(J,II)=A(J,I) 117 CONTINUE DO 118 II=1,M I=IPT(II) A(II,II)=D(I) 118 CONTINUE C RETURN END C C____________________________________________________________________________ C SUBROUTINE SORTC(N,M) IMPLICIT INTEGER*2 (I-N) PARAMETER (JMAX=250,MAXDIM=JMAX/2+1,NLINES=10000) C COMMON /SORTCC/WK,IPT INTEGER*2 IPT(NLINES) REAL*8 WK(NLINES),EE C C ... This routine sorts the quantities in vector WK in ascending order C of magnitude and also accordingly rearranges vector IPT of pointers C to original positions of