C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C QSTARK - PROGRAM TO FIT THE STARK SHIFTS FOR A LINEAR/SYMMETRIC/ASYMMETRIC C ROTOR WITH UP TO ONE QUADRUPOLAR NUCLEUS + SPIN-ROTATION C C (this is modified from a two quadrupole I,F coupling program, C - the two quadrupole mechanisms have been made dormant but C will hopefully be reactivated in the future) C C Refs: C C (H2O)2...HCl - Z.Kisiel, J.Kosarzewski, B.A.Pietrewicz, L.Pszczolkowski, C Chem.Phys.Lett. 325, 523-530 (2000) C H2O...HCl - Z.Kisiel, B.A.Pietrewicz, P.W.Fowler, A.C.Legon, E.Steiner, C J.Phys.Chem.A 104,6970-6978 (2000) C C C Ver 10.I.2023 ------- Zbigniew KISIEL ------- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, 02-668 Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c C Modification history: C C 29.07.99: first working version developed from Q2FIT C 20.09.99: modification to fit frequency differences C 30.09.99: frequency contributions made operational C 19.03.00: some modifications to output + field correction constant k C 30.03.00: changes to output/input C 23.07.00: centrifugal distortion on chi C 14.08.00: improvements to eigenvalue assignment C 27.09.00: elimination of longstanding bug in LSQ which affected deviations C on constants when there was a large number of unfitted lines C 18.11.00: predictive mode and diagnostic plots C 26.11.00: gle output C 8.12.00: vertical gle output+more gle mods. C 12.02.01: incremental mods including input checks C 24.08.01: mods for high F and high MDIM operation C 5.12.01: inclusion of the Chi.d term C 9.12.02: proper treatment of mu.c as an explicit imaginary term C 6.09.05: correction of bug in 'frequency differences' mode and mods C 30.01.07: modification for a weighted fit C 30.07.07: dynamic dimensioning of the largest arrays C 18.09.07: various mods and elimination of bug affecting deltaC C 10.10.07; option of relative or absolute frequencies in gle output C 24.08.09: completion of treatment of chi_bc and chi_ac C 22.12.09: modified gle and contributions output C 4.10.17: porting of CONFIX from PIFORM for clearer printout of fixed parameters C 19.05.20: eliminating obsolete Fortran constructs and debugging c 28.07.20: tweaks to gle output c 13.08.20: default output of EXP.DAT when in the fitting mode c 26.08.22: worst fitted transitions + more internal commenting c 14.11.22: write EXPTPLOT.DAT also when in fitting mode+debugging c 10.01.23: less truncated file names for gle output C C_____________________________________________________________________________ C C C The algorithm is as follows: C C -------- Main fitting loop C | C | ------ Loop over distinct field values C | | C | | ---- Loop over distinct M.F values C | | | C | | | 1. Set up and diagonalize H (real or imaginary, as required) C | | | 2. Fish out the appropriate eigenvalues and set up C | | | derivatives C | | |---------------------- C | -------------------------------- C | Evaluate a new set of constants C ---------------------------------------- C C Recommended use is to: C C 1/ First run the program to fit the electrode spacing D from calibration C measurements C 2/ Then use D from above to fit dipole moment components C C Make sure that the NQOVER parameter (see below) is not too small - the C eigenmatrix is infinite and NQOVER controls truncation. Special C consideration for higher symmetries: C C LINEAR TOP: set A to some suitably large value C set B to B C set C to zero C use asymmetric top labelling for a-type K_a=0 C transitions and fit the mu_a dipole component C C SYMMETRIC TOP: set A to the value of the axial constant (C for oblate top) C set B to B C set C to zero C C In addition, for symmetric tops: C C - Use asymmetric top level labelling, and for K>0 levels of C symmetric tops the lower energy level will be that with lower C Tau. Since Stark splitting decreases with J this means that C Stark components shifted to high frequency will tend to be C associated with lower values of Tau and vice versa. C - The fitted (axial) dipole moment is to be specified as mu_a, C irrespective of whether it is in reality mu_a or mu_c C C C PROGRAM FEATURES: C C - All Delta_M transition types (0, -1, +1) can be fitted simultaneously C C - Enhanced determination of spectroscopic parameters from local C perturbations is possible (up to the limit of applicability of the C eigenvalue identification scheme - problems will appear near level C crossings) C C - Calculation for tops without quadrupolar nuclei can be made C and results of various orders can be fitted simultaneously C C - possibility of a weighted fit C C - Calculations for all types of tops, asymmetric, symmetric, and linear, C with or without a quadrupolar atom have been tested C C - Prediction of Stark lobe behaviour over a specified voltage range, C for both linear and quadratic dependence C C - ASCII diagnostic plots of Stark lobe behaviour C C - gle plots of Stark lobe behaviour C C - Field correction constant k is available for empirical correction of C departure of field from the applied voltage, through polarization or C other effects: the actual field is given by V(1+kV)/d. k may be C required if inclusion of successively higher field measurements is C found to systematically shift the dipole moment (but note that for C high fields polarizability effects have to be accounted for) C C C KNOWN BUGS/DEFICIENCIES: C C The program calculates correct energies but is known to run into labelling C problems when the off-diagonals in the H matrix become sufficiently large. C Known instances of such behaviour are: C C - first order Stark effects in symmetric tops C - highly perturbing states in asymmetric tops C - high field calculations when mu_c is non-zero and QSTARK switches to C the complex H matrix form C C Extensive dump output, as controlled by the IDUMP parameter, allows C more detailed insight into such behaviour. C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C COUPLING SCHEME: F = J + I C C FOR MATRIX ELEMENTS SEE: C C H.E H.P.Benz,A.Bauder,Hs.H.Gunthard - J.Mol.Spectrosc. 21,156-164(1966) C C H.Q Keenan,Wozniak,Flygare - J.Chem.Phys. 75,631(1981) C Centrifugal distortion on chi is defined using C Chi.aa = chi.aa+ chi.j J(J+1) + chi.k K**2 etc. C The above and the Chi.d constant originate from C J.T.Hougen - J.Chem.Phys. 57,4207(1972) C Chi.d is implemented as a simple symmetric top-type diagonal C contribution as used (and reproducing the results) in: C Burie,Boucher,Demaison,Dubrulle - Molec.Phys. 32,289(1976) C Fowler,Legon,Thumwood,Waclawik - Coord.Chem.Rev. 197,231(2000) C C H.SR Read,Flygare - J.Chem.Phys. 76,2238(1982) C Note that the SR terms are defined for H.sr=-I.M.J and are C of opposite sign to the analogous constants in Pickett's programs C which are for H.sr=+I.M.J C C FOR WEIGHTED LEAST SQUARES FITTING SEE: C C D.L.Albritton,R.N.Schmeltekopf,R.N.Zare - "An introduction to the C least-squares fitting of spectroscopic data", pp.1-66 in "Molecular C Spectroscopy: Modern Research", Vol.2, Academic Press (1976) C C C QSTARK has been developed from program Q2FIT. Many of the original C comments pertaining to Q2FIT are still in place and some may no C longer be applicable. The comments immediately following are old c comments from Q2FIT: C C C This program is a fitting version of modified program READ written C by G.W.Read ca 1981, containing angular momentum routines written C by M.P.Keenan ca 1980. C Some enhancements to elements are by Claudio di Esposti. C C NOTE: Q2FIT has been exhaustively tested for the single nucleus C case including X.ab. The calculation for two nuclei has been tested C with Ar...ClCN. Spin rotation has only been tested in the C symmetric limit, and spin-spin has not yet been tested at all. C C ERRATA to papers: C C H.Q: C 1/ In eq.11 lower line second term (the 3j symbol following the C 6j symbol) should have -q instead of -2 in its lower line, so C that it should read J 2 J' C -K -q K' C 2/ The last two spherical tensor components in eq.12 should have C signs -+ and + resp. (instead of +- and +-). C 3/ The off-diagonal quadrupole constants are actually C X.ab+-iX.ac and X.bb-X.cc+-2iX.bc C (QSTARK currently uses only the single constant X.ab) C H.SR: C 1/ in eq.A3 the closing square bracket after delta(J'J)delta(K'K) C is missing. C C_____________________________________________________________________________ c C C DATA FILE IS TO BE REEDITED FROM A PREVIOUS EXAMPLE which should be C self-explanatory. C c///////////////////////////////////////////////////////////////////////////// C HEADER: Data set begins with a header containing c///////////////////////////////////////////////////////////////////////////// C C COMENT - up to 78 character alphanumeric comment C dummy line - obligatory but not used by the program C ITYPF - defines the type of fit: C =1 fit directly to transition frequencies C =2 fit to frequency differences from appropriate field free C transitions C I1,I2 - spins of the two quadrupolar nuclei (in units of 1/2), both C can be set to zero for spectrum without hyperfine structure C (in this version I2 is always zero) C NIT - mode selection, translated internally into NPREDS, C >0 fitting mode, number of iterations (NPREDS=0) C =0 predictive mode (NPREDS=1) C =-1 and -2 predictive mode with landscape gle output C -1 plot of frequency.vs.voltage (NPREDS=2) C -2 plot of frequency.vs.field (NPREDS=3) C =-11 and -12 as -1 and -2 except for portrait orientation C of the gle plot C NOTE: C a/ when STARK predictive mode gle output is selected it is C possible to place measured data points on top of the generated C curves. In order to do this it is necessary to put C a simple two column list of voltage,frequency C pairs in file EXP.DAT. This file can be generated automatically C during a predictive run, or can be edited/reedited with text C editor. When QSTARK prediction is run this file is C converted to EXPTPLOT.DAT, which will then be used in C the generation of PDF/PostScript/JPG plot. C IMPORTANT: C If this is to work smoothly for a plot of C frequency differences then it is necessary C for EXP.DAT to contain points in such a way that C they are in blocks for Stark components, each C starting with zero field frequency for this C component. C b/ When QSTARK is run in the fitting mode then it will C automatically generate the EXP.DAT file for possible C subsequent use in the predictive mode. C The EXPTPLOT.DAT file is also written at the end of C a fitting run, but with the limitation that the E and E**2 C columns are always voltages (i.e. cell constant d=1). C This file can be used directly with a previously generated C gle diagram without repeating the predictive run. C This allows a quick check of the positions of newly measured C points, though at the cost of using a perhaps outdated C prediction. C C NQOVER - the excess in F and J, over the values for the energy level of C interest, to which the Hamiltonian is to be set up C (in general NQOVER=0 is unsatisfactory, for most cases C NQOVER=1 yields results as good as NQOVER>1, and for Br it C is best to use NQOVER=2) C IDUMP - allows dumping of additional information on the internal C working of the program to a file called DUMP, and sometimes C to screen C =0 no dumps C =1 all raw eigenvalues as obtained out of diagonalization + C results of subsorting within J blocks C =2 as above + mapping of raw eigenvalues onto J blocks C + screen output of various indices and dimensions C =3 as above + full hamiltonian matrix C ICONTR - allows assessment of contributions to frequencies C =0 none, standard output C =1 contributions of fitted constants to frequencies, these C are evaluated from gradients so that they will be C subject to nonlinearities C =2 field free frequencies, this option is really a leftover C from quadrupole free option of Q2FIT but is useful in C assessing the nonlinearity in stark shift curves since the C calculation is F.free=F.obs-sum(derivative*constant) C NT - the number of transitions in the data set. In the case of C a standard fitting deck this can be set to some large value C and the program will count the input lines automatically, C providing the string 'end' is placed in a line directly after C the last data line C For predictions NT has to be declared exactly. C dummy line - obligatory but not used by the program C c c///////////////////////////////////////////////////////////////////////////// C CONSTANTS: NCONST lines declaring the values C of spectroscopic constants and their fitting parameters c///////////////////////////////////////////////////////////////////////////// C C--------line declaring a constant: C C DUMMY IA(1) CON(1) :FORMAT (A14,I2,F30) C C DUMMY - 14 character alphanumerical variable allowing for use C of constant names in the data file (this is not used C by the program, which generates the names internally) C IA(i) - fitting parameters for the spectroscopic constants C (=1 constant is to be fitted, =0 fixed) C CON(i) - spectroscopic constants C c The constants (Watson's A-reduced Hamiltonian and Flygare's C H_sr=-I.M.J convention for signs of M): c c A,B,C - rotational c DJ,DJK,DK,deltaJ,deltaK - quartic c.d. c HJ,HJK,HKJ,HK,hj,hjk,hkj - sextic c.d. c LKKJ - octic c.d. c Chi_aa,Chi_bb-Chi_cc - determinable diagonal chis c Chi_ab - the most common off-diagonal chi c Chi^J_aa,Chi^K_aa,Chi^J_b-c - centrifugal distortion of diagonal chis c M_aa,M_bb-M_cc,Trace(M) - spin-rotation c Chi_d - Hougen's chid constant c Chi_bc,Chi_ac - remaining two off-diagonal constants c nu,nu,nu - three constants that are not yet used c mu_a,mu_b,mu_c - dipole moment components c d - cell constant c k - empirical correction for nonlinearity of c cell calibration with field c c c///////////////////////////////////////////////////////////////////////////// c TRANSITIONS: NT lines declaring c 1/ either Stark and/or field free measurements (fitting mode) c 2/ or Stark lobes (predictive mode) c///////////////////////////////////////////////////////////////////////////// C C 1/-----line declaring a Stark measurement (F and MF are to be in units of 1/2): c C J' K_a' K_c' J K_a K_c F' M.F' F M.F Frequency VOLTS LINFIT ERROR C C which is read in free format. C C FREQUENCY - frequency of transition i, which is placed in TFREQ(i,1) C VOLTS(i) - voltages at which the transitions have been observed C (in the same order as frequencies in TFREQ), the voltages are C scaled into fields using the electrode spacing, D, which is C the last of the declared constants C LINFIT(i) - pointer to whether line i is to be fitted (=1) or omitted C from the fit, ie predicted only (=0) C ERROR(i) - error in measured transition frequency that will be used for C weighting purposes (wt=1/error^2). C This value is optional, if it is not provided then the error C will have the value of zero and the weight will be set to unity. C C --> The number of lines can be counted by the program - to do this C set NT to some large value and place string 'end' in a line C directly after the last data line C (ie. just the three characters e,n,d beginning in the C first column) C C --> Empty lines can be placed between declaration lines, and will C be echoed as such in the same places in the printout C C --> Lines between transition declaration lines beginning with C the '!' or '%' character are treated as comment lines - those C beginning with '!' will be echoed in the output C C C 2/-----line declaring a Stark component for prediction as a function of C applied voltage: C C J' K_a' K_c' J K_a K_c F' M.F' F M.F C Frequency VOLMIN VOLMAX ITYPRE NVSTEP FLOWER FUPPER C C which is read in free format. C C VOLMIN - minimum value of voltage for use in predicting voltage C behaviour of a Stark lobe C VOLMAX - maximum value of voltage for use in predicting voltage C behaviour of a Stark lobe C ITYPRE =1 divide the voltage into linear increments C =2 divide the voltage to obtain linear increments in V**2 C NVSTEP - number of voltage increments up to VOLMAX C FLOWER - lower limit for Stark shift in diagnostic plot C FUPPER - upper limit for Stark shift in diagnostic plot C If both FLOWER and FUPPER are set to zero (put in three commas) C then the vertical axis of the plot is autoscaled. C Only the FLOWER,FUPPER parameters of the last data line C matter but, owing to free format input, at least the C appropriate separators such as ,,, have to be placed in C the data line. C For absolute frequency plot FLOWER,FUPPER are applied to C the mean zero field frequency of specified components. C C --> Declarations of Stark components do not allow insertion of C comments between the lines C --> Automatic counting of the number of declared Stark lobes is C not implemented so NT has to be declared exactly C C C QSTARK asks for input and output file names, but these can be C put in a separate file, say QS.IN, for pipelined operation of the C type QSQS.OUT &. C C_____________________________________________________________________________ C C C FH,FL - highest and lowest values of the quantum number F (in units C of 1/2) C MFH,MFL - highest and lowest values of the quantum number MF (in C units of 1/2) C NCON - number of constants to be fitted C ICFIT(i) - pointers to constants to be fitted (ie. only the first NCON C elements are set C PO(i) - names of constants used in the printout C POSHRT(i) - short names of constants for correlation coefficients and C contributions parts of printout C X(i) - changes in values of constants after a least squares iteration C SD(i) - standard deviations on fitted constants C C LABO(i,j,k) - quantum numbers for transitions: C i is the number of the transition C j defines the level involved in the transition and is 1 for the C lower and 2 for the upper C k is a quantum number index for a level in the transition: C =1 I C =2 J C =3 Tau+J+1 which runs from 1 to 2J+1 for a given J, and counts C the position of the energy level in sorted levels in a C given J block C =4 F C =5 M.F C TFREQ(i,j) - frequencies of transitions for i from 1 to NT where: C (i,1) observed frequencies C (i,2) calculated frequencies C (I,3) obs-calc C IWFIT - specifies whether the fit is to be equally weighted C or weighted on the basis of declared experimental errors, C which are optionally placed after the fitting parameter of C each declared transition. IWFIT is set internally on C the basis of whether errors have been declared or not: C IWFIT=0 equally weighted fit (no errors declared) C IWFIT>0 weighted fit (errors have been found in the input, in C which case each line should carry such an error) C ERROR(i) - error in frequency i C WT(i) - weight for fitting frequency i = 1/ERROR(i)^2 C REFERR(i) - copy of error in frequency i as read in from the data file C for use in evaluating the errors in frequency differences c C INFOE(i,j) - information on energy levels to be calculated where: C (i,1) points to the transition the level is involved in and C is negative for the lower, and positive for the upper level C in the transition C (i,2) is the value of M.F for the level (in units of 1/2) C VOLINF(i) - voltages for energy levels defined in INFOE, in the same order C INFTMP(i,j) - temporary copy of INFOE used for sorting purposes C IFREE(i) - pointers to field free transitions associated with Stark C component i C INDEX(i) - index built up from J and Tau quantum numbers of both levels C in transition i C INDEXF(i) - index built up from F quantum numbers of both levels in C transition i: the values of INDEX and INDEXF provide C unique identification of the transition C C COMBLK - string variable holding all comments 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 C C_____________________________________________________________________________ C C Compilation: C C Static allocation of local variables is crucial for this program. This is C no longer a default in contemporary Fortran compilers, so that you have to C enforce it with appropriate keywords: C C Compaq Visual Fortran: -static C Intel Visual Fortran: -Qsave C C///////////////////////////////////////////////////////////////////////////// C C C Declaration of all large arrays for global use by all program units C MODULE globdat integer(4) mdim,msize C C H1(i,j) - Hamiltonian matrix C T(i,j) - eigenvectors from diagonalisation of H1 C NELINF(i,j) - information on eigenvalues of H required by the spectral data: C (i,1) is the number of transition the level belongs to (-ve for C the lower level) C (i,2) is the number of this eigenvalue among the eigenvalues of H C LINSYM - top symmetry flag, =0 for asymmetric, =1 for linear and C symmetric molecules C REAL(8) H1 REAL(4) T INTEGER(4) NELINF,LINSYM c ALLOCATABLE :: H1(:,:),T(:,:),NELINF(:,:) C END MODULE globdat C C///////////////////////////////////////////////////////////////////////////// C C C M A I N P R O G R A M c use globdat c IMPLICIT REAL(8) (A-E,G-H,O-Z), INTEGER(F) PARAMETER (NCONST=36,NLIN=2000,NLIN2=2*NLIN,msort=10000, * nerd=2,nvpix=52,nhpix=79) PARAMETER (maxcbl=10000,maxcom=200,maxspa=200) c COMMON /BLKEN/TFREQ,DER,LABO,INFOE,volinf,ICFIT,NELEV COMMON /BLKFIT/WT,LINFIT,IFREE(nlin),CORELC(NCONST,NCONST) COMMON /sortcc/wk,IPT common /stark/volts(nlin) common /cmnt/coment,descdat DIMENSION CON(NCONST),IA(NCONST),volinf(nlin2) DIMENSION LABO(NLIN,2,5),DER(NLIN,NCONST),LINFIT(NLIN) DIMENSION ICFIT(NCONST),TFREQ(NLIN,3),INFOE(NLIN2,2), * inftmp(nlin2,3),ERROR(nlin),WT(nlin),REFERR(nlin) DIMENSION SD(NCONST),X(NCONST),INQ(0:NLIN,10) INTEGER(4) index(0:nlin),indexf(0:nlin),UNITS(NCONST) dimension wk(msort),ipt(msort),ispace(0:maxspa),icomm(maxcom), * icomml(maxcom) CHARACTER COMENT*78,comcop*80,line*80,errmes*30,comblk*(maxcbl), * DESCDAT*59(200) CHARACTER FILNAM*30,filout*30,linout*140,worstlin(6)*140, * HEADER1*140,HEADER2*140 CHARACTER DUMMY*14,PO(NCONST)*12,POSHRT(NCONST)*6 CHARACTER CONVAL*25,ERVAL*25,plotbf(nhpix),cbuff*50 CHARACTER OUTCON(NCONST)*39,clabel*15 REAL(8) fmin,fmax,flow,fhigh,fincr,flower,fupper,fmark,fcurm, * fzmin,fzmax,fzmean,worstomc(6),fzero CHARACTER(7) UNIT,NAMUN(8) CHARACTER(96) FORMAL,F1,F2 EQUIVALENCE (cbuff,plotbf(1)),(clabel,plotbf(54)), * (comcop,plotbf(1)) 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 ','D ','cm '/ DATA UNITS/1,1,1, 2,2,2,2,2, 3,3,3,3,3,3,3, 4, 1,1,1, * 2,2,2, 1,1,1, 2,1,1, 2,2,2, 7,7,7, 8,8/ C DATA PO/' A ',' B ',' C ', ! 1- 3 1 ' DJ ',' DJK ',' DK ', ! 4- 6 1 ' dJ ',' dK ',' HJ ', ! 7- 9 1 ' HJK ',' HKJ ',' HK ', ! 10-12 1 ' hJ ',' hJK ',' hK ', ! 13-15 1 ' LKKJ ', ! 16 1 'I1: X.aa ','I1: Xbb-Xcc ','I1: Xab ', ! 17-19 1 'I1: XJ.aa ','I1: XK.aa ','I1: XJ.b-c ', ! 20-22 1 'I1: Maa ','I1: Mbb-Mcc ','I1: Trace(M)', ! 23-25 1 'I1: X.d ','I1: Xbc ','I1: Xac ', ! 26-28 1 ' ',' ',' ', ! 29-31 1 ' Mu.a ',' Mu.b ',' Mu.c ', ! 32-34 1 ' d ',' k '/ ! 35-36 c DATA POSHRT/'A ','B ','C ', ! 1-3 1 'DJ ','DJK ','DK ','dJ ','dK ', ! 4-8 1 'HJ ','HJK ','HKJ ','HK ','hJ ','hJK ','hK ',! 9-15 1 'LKKJ ', ! 16 1 '1:Xa ','1:Xb-c','1:Xab ', ! 17-19 1 '1:XJ.a','1:XK.a','1:XJbc', ! 20-22 1 '1:Ma ','1:Mb-c','1:Tr ', ! 23-25 1 '1:Xd ','1:Xbc ','1:Xac ',' ',' ',' ', ! 26-31 1 'Mu.a ','Mu.b ','Mu.c ',' d ',' k '/ ! 32-36 c F1='(I3,1H., 1x,3I3,1x,3I3,2X,2(I3),1x,2(I3),F10.1,F14.4,3x'// * ',F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' F2='(66x '// * ',F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0,F9.0)' C C C...Input C Write(*,'(1x//)') WRITE(*,200) 200 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | QSTARK - Fully diagonalizing fit of Stark shifts in a', * ' rotor',T79,'|'/ * ' | with zero or one quadrupolar nuclei', * T79,'|'/ * ' |',76(1H_),'|'/' version 10.I.2023', * T64,'Zbigniew KISIEL'//) c 202 WRITE(*,201,advance='NO')'Name of data file: ' 201 FORMAT(1X/1X,A) READ(*,'(A)',ERR=202)FILNAM OPEN(2,FILE=FILNAM,ERR=202,status='old') READ(2,'(A)')COMENT read(2,'(a)')cdummy ! skip separating line READ(2,260,err=261)Itypf if(itypf.eq.1.or.itypf.eq.2)goto 262 c 261 write(*,'(1x,a,i5/)')'Incorrect type of calculation flag: ', * itypf stop c 262 READ(2,260,err=263)I1 read(2,260,err=263)I2 I2=0 goto 264 c 263 write(*,'(1x,a/)')'Problem with nuclear spin: ' stop c 264 READ(2,260,err=265)NIT goto 266 c 265 write(*,'(1x,a/)')'Problem with number of iterations' stop c 266 if(NIT.gt. 0)then npreds=0 else if(nit.lt.-10)then nit=nit+10 iveho=1 else iveho=0 endif npreds=iabs(nit)+1 endif c read(2,260,err=269)NQOVER read(2,260,err=269)IDUMP read(2,260,err=269)ICONTR goto 270 c 269 write(*,'(1x,a/)')'Problem reading NQOVER,IDUMP or ICONTR' stop c 270 READ(2,260,err=267)NT 260 format(38x,i3) goto 268 c 267 write(*,'(1x,a/)')'Problem with number of transitions' stop c 268 read(2,'(a)')cdummy ! skip separating line C C...Read parameters C NCON=0 DO 350 I=1,NCONST read(2,'(a)')line READ(line,370,err=390)DUMMY,IA(I),CON(I) 370 FORMAT(A14,I2,F30.0) IF(IA(I).EQ.0)GOTO 350 NCON=NCON+1 ICFIT(NCON)=I 350 CONTINUE goto 391 c c...patch to avoid a crash on older data without the k parameter c 390 if(i.eq.36)then ia(i)=0 con(i)=0.0d0 backspace(2) else write(*,'(//1x,a/1x,a//)') * 'Problem reading a parameter from input line:', * line(1:len_trim(line)) stop endif c 391 linsym=0 if(con(3).eq.0.d0)linsym=1 c c...do not fit any constants in predictive mode c if(npreds.ge.1)then do 42 i=1,nconst ia(i)=0 42 continue ncon=0 endif C C...order of quantum numbers: C 1 2 3 4 5 6 7 8 9 10 C C LABO: I J Tau+J+1 F M.F C INQ: J' K_a' K_c' J K_a K_c F' M.F' F M.F C c c...Read lines for fitting c if(npreds.eq.0)then nzero=0 notzer=0 NLFIT=0 nspace=0 ncomm=1 icomm(ncomm)=0 IWFIT=0 DO 2 I=1,NT errmes='unreadable line or end of file' 23 read(2,'(a)',err=351,end=351)line 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)=i icomm(ncomm)=icomm(ncomm-1)+lentrm(line) comblk(icomm(ncomm-1)+1:icomm(ncomm))=line(1:lentrm(line)) goto 23 endif c if(line(1:3).eq.'end')then NT=i-1 goto 22 endif c if(lentrm(line).lt.20)then nspace=nspace+1 if(nspace.lt.maxspa)ispace(nspace)=i if(nspace.gt.1.and.i.eq.ispace(nspace-1))nspace=nspace-1 goto 23 endif c READ(LINE,*,ERR=500,end=500)(INQ(I,J),J=1,10),TFREQ(I,1), * volts(i),LINFIT(I),error(i) c if(error(i).le.0.d0)then error(i)=0.d0 wt(i)=1.d0 else wt(i)=1.d0/error(i)**2 IWFIT=IWFIT+1 endif referr(i)=error(i) goto 501 c 500 READ(LINE,*,ERR=351,end=351)(INQ(I,J),J=1,10),TFREQ(I,1), * volts(i),LINFIT(I) error(i)=0.d0 wt(i)=1.d0 c 501 errmes='illegal quantum numbers' isum=inq(i,1)-inq(i,2)-inq(i,3) if(isum.ne.0.and.isum.ne.-1)goto 351 isum=inq(i,4)-inq(i,5)-inq(i,6) if(isum.ne.0.and.isum.ne.-1)goto 351 if(inq(i,7).gt.2*inq(i,1)+I1)goto 351 if(inq(i,9).gt.2*inq(i,4)+I1)goto 351 if(inq(i,7).lt.iabs(2*inq(i,1)-I1))goto 351 if(inq(i,9).lt.iabs(2*inq(i,4)-I1))goto 351 if(inq(i,8).gt.inq(i,7))goto 351 if(inq(i,10).gt.inq(i,9))goto 351 c IF(LINFIT(I).EQ.1)then NLFIT=NLFIT+1 if(volts(i).eq.0.d0)then nzero=nzero+1 else notzer=notzer+1 endif endif LABO(I,1,1)=I1 LABO(I,1,2)=INQ(I,4) LABO(I,1,3)=INQ(I,5)-INQ(I,6)+INQ(I,4)+1 LABO(I,1,4)=INQ(I,9) LABO(I,1,5)=INQ(I,10) LABO(I,2,1)=I1 LABO(I,2,2)=INQ(I,1) LABO(I,2,3)=INQ(I,2)-INQ(I,3)+INQ(I,1)+1 LABO(I,2,4)=INQ(I,7) LABO(I,2,5)=INQ(I,8) index(i)=labo(i,1,2)*1000000+labo(i,1,3)*10000+ * labo(i,2,2)*100 +labo(i,2,3) indexf(i)=labo(i,1,4)*100+labo(i,2,4) 2 CONTINUE 22 GOTO 352 endif c c...read lines for prediction c if(npreds.ge.1)then c nreadl=0 nlfit=0 27 nlfit=nlfit+1 i=nreadl+1 errmes='unreadable line or end of file' read(2,'(a)',err=351,end=351)line backspace(2) READ(2,*,ERR=351,end=351)(INQ(nlfit,J),J=1,10),TFREQ(nlfit,1), * volmin,volmax,itypre,nvstep,flower,fupper c errmes='error in prediction parameters' if(nvstep.le.0)goto 351 if(volmax.le.volmin)goto 351 if(fupper.lt.flower)goto 351 if(itypre.lt.1.or.itypre.gt.2)goto 351 c errmes='illegal quantum numbers' isum=inq(nlfit,1)-inq(nlfit,2)-inq(nlfit,3) if(isum.ne.0.and.isum.ne.-1)goto 351 isum=inq(nlfit,4)-inq(nlfit,5)-inq(nlfit,6) if(isum.ne.0.and.isum.ne.-1)goto 351 if(inq(nlfit,7).gt.2*inq(nlfit,1)+I1)goto 351 if(inq(nlfit,9).gt.2*inq(nlfit,4)+I1)goto 351 if(inq(i,7).lt.iabs(2*inq(i,1)-I1))goto 351 if(inq(i,9).lt.iabs(2*inq(i,4)-I1))goto 351 if(inq(nlfit,8).gt.inq(nlfit,7))goto 351 if(inq(nlfit,10).gt.inq(nlfit,9))goto 351 c nreadl=nreadl+1 nlfit=nlfit-1 do 25 ii=1,nvstep+1 nlfit=nlfit+1 if(itypre.eq.1)then volts(nlfit)=volmin+(ii-1)*(volmax-volmin)/real(nvstep) else volts(nlfit)= * dsqrt( volmin**2+(ii-1)*(volmax**2-volmin**2)/real(nvstep) ) endif linfit(nlfit)=0 c i=nlfit do 26 j=1,10 inq(i,j)=inq(i-(ii-1),j) 26 continue tfreq(i,1)=tfreq(i-(ii-1),1) LABO(I,1,1)=I1 LABO(I,1,2)=INQ(I,4) LABO(I,1,3)=INQ(I,5)-INQ(I,6)+INQ(I,4)+1 LABO(I,1,4)=INQ(I,9) LABO(I,1,5)=INQ(I,10) LABO(I,2,1)=I1 LABO(I,2,2)=INQ(I,1) LABO(I,2,3)=INQ(I,2)-INQ(I,3)+INQ(I,1)+1 LABO(I,2,4)=INQ(I,7) LABO(I,2,5)=INQ(I,8) index(i)=labo(i,1,2)*1000000+labo(i,1,3)*10000+ * labo(i,2,2)*100 +labo(i,2,3) indexf(i)=labo(i,1,4)*100+labo(i,2,4) if(nlfit.ge.nlin)goto 28 25 continue c if(nreadl.lt.NT)goto 27 c 28 NT=nlfit nlfit=0 NIT=1 goto 352 endif c 351 WRITE(*,353)i,errmes,line 353 format(1X/' INPUT ERROR ON LINE',I5,' ----> ',a/ * ' check:'/1x,a/) STOP c 352 CLOSE(2) C if(npreds.eq.0)then WRITE(*,257)NT,NLFIT,nzero,notzer,NCON 257 FORMAT(1X/1X,I5,' Lines read in'/1X,I5,' Lines fitted --->', * i5,' field free and',i5,' at non-zero field' * /1X,I5,' Constants fitted') else WRITE(*,1257)NT,nlin 1257 format(1x/1x,I5,' Lines have been generated (maximum is',i5,')') endif NDEGF=NLFIT-NCON IF(NDEGF.LT.1)NDEGF=1 c c c...set up pointers for fitting frequency differences C if(itypf.eq.2)then do 403 I=1,NT ifree(i)=0 403 continue c do 400 I=1,NT c if(volts(i).eq.0.0)then ifree(i)=i if(tfreq(i,1).eq.0.0)then write(*,'(1x/'' Field free transition:'',i4, * '' has zero frequency''// * '' INPUT ERROR: Please correct the data set''//)')i stop endif else goto 400 endif c do 401 ii=1,nt if(ii.eq.i)goto 401 if(( index(i).eq. index(ii)).and. * (indexf(i).eq.indexf(ii)) )then ifree(ii)=i endif 401 continue c 400 continue c ii=0 do 402 I=1,NT if(ifree(i).eq.0)then if(ii.eq.0)write(*,'(1x//)') write(*,'('' Transition:'',i4,'' at '',f12.5, * '' MHz and'',f8.1,'' Volts has no field free value'')') * i,tfreq(i,1),volts(i) ii=ii+1 endif 402 continue if(ii.ne.0)then write(*, * '(1x//'' INPUT ERROR: Please correct the data set''//)') stop endif endif c C...Initialize output C 203 WRITE(*,201,advance='NO')'Name of output file: ' READ(*,'(A)',ERR=203)FILOUT OPEN(3,FILE=FILOUT,ERR=203,STATUS='UNKNOWN') write(3,200) WRITE(3,1)COMENT,nqover 1 FORMAT(1X,78(1H-)/1X,A/1X,78(1H-)// * ' Calculation for J and F exceeding the value in data by', * ' at least',I2) if(npreds.eq.0)then if(itypf.eq.2)then write(3,410)'FREQUENCY DIFFERENCES' else write(3,410)'TRANSITION FREQUENCIES' endif WRITE(3,257)NT,NLFIT,nzero,notzer,NCON else write(3,41) write(3,1257)NT,nlin endif 410 format(' The fit will be made to ',a) 41 format(' Prediction of STARK LOBE VOLTAGE DEPENDENCE') C c C...Extract the energy levels to be calculated and then sort them: c first - in ascending order of voltage C second - within each voltage value in ascending order of M.F C NELEV=0 DO 230 I=1,NT DO 231 J=1,2 NELEV=NELEV+1 INFOE(NELEV,1)=-I IF(J.EQ.2)INFOE(NELEV,1)=I INFOE(NELEV,2)=LABO(I,J,5) volinf(nelev)=volts(i) 231 CONTINUE 230 CONTINUE C c...sort energy levels in INFOE according to increasing voltage c if(nelev.gt.msort)then write(*,'(1x// * '' ERROR: the current MSORT='',i5,'' is too small''/ * '' and needs to be increased''//)') * msort stop endif c do 300 j=1,nelev wk(j)=volinf(j) ipt(j)=j inftmp(j,1)=infoe(j,1) inftmp(J,2)=INFOE(J,2) 300 continue call sortc(1,nelev) ! <---- SORTC do 301 j=1,nelev jj=ipt(j) infoe(j,1)=inftmp(jj,1) infoe(j,2)=inftmp(jj,2) volinf(j)=wk(j) 301 continue c c...subsort levels for each voltage in order of M.F c jstart=1 ndistv=0 304 vstart=volinf(jstart) do 302 j=jstart,nelev if(volinf(j).ne.vstart)then jfin=j-1 goto 303 endif 302 continue c jfin=nelev 303 ndistv=ndistv+1 if(jfin-jstart.eq.0)then jstart=jstart+1 if(jstart.ge.nelev)goto 306 goto 304 endif do 305 j=jstart,jfin ipt(j)=j wk(j)=infoe(j,2) inftmp(j,1)=infoe(j,1) inftmp(j,2)=infoe(j,2) 305 continue call sortc(jstart,jfin) ! <---- SORTC do 310 j=jstart,jfin jj=ipt(j) infoe(j,1)=inftmp(jj,1) infoe(j,2)=inftmp(jj,2) 310 continue jstart=jfin+1 if(jstart.lt.nelev)goto 304 306 if(idump.gt.1)then write(*,307)nelev,ndistv write(3,307)nelev,ndistv 307 format(/i6,' energy levels'/i6,' distinct voltages'// * ' ntrans M.F V I J t+J+1 F M.F'/) do 308 j=1,nelev labpt=2 if(infoe(j,1).lt.0)labpt=1 write(*,309)infoe(j,1),infoe(j,2),volinf(j), * (labo(iabs(infoe(j,1)),labpt,jj),jj=1,5) write(3,309)infoe(j,1),infoe(j,2),volinf(j), * (labo(iabs(infoe(j,1)),labpt,jj),jj=1,5) 308 continue endif 309 format(2i5,f10.2,5x,5i6) c c if(idump.ge.1)open(4,file='dump',status='unknown') ! dump C C C------------------ MAIN FITTING LOOP -------------------- C C N=0 DO 15 LP=1,NIT c do 465 n=1,5 worstomc(n)=0.d0 465 continue C if(idump.ne.0)then ! dump do 1300 i=1,NCONST ! dump if(con(i).eq.0.d0)goto 1300 ! dump write(4,1301)po(i),con(i) ! dump 1301 format (1X,A,' =',F20.9) ! dump 1300 continue ! dump endif ! dump c DO 3000 J=1,NT ! Fortran2018 TFREQ(J,2)=0.D0 DO 3 M=1,NCON DER(J,M)=0.0D0 3 CONTINUE 3000 CONTINUE ! Fortran2018 C C...set up H, diagonalise, and calculate frequencies C WRITE(*,210)'CALL TO READ - Iteration',LP 210 FORMAT(1X/1X,A,I3) C CALL READ(CON,NCON,I1,I2,IDUMP,nt,nqover) ! <---- READ C C...write output header C write(3,1256) 1256 format(1X/' TRANSITIONS (F and MF in units of 1/2):'//) c if(itypf.eq.1.and.iwfit.eq.0)then ! freq fit,unweighted WRITE(3,255) WRITE(3,1255) write(header1,255) write(header2,1255) endif 255 FORMAT( * 6X,'J K K <- J K K F MF<- F MF',4X,'Volts',8x, * 'Obs',10X,'Obs-Calc',' Calc ',6x, * 'Field Field**2') 1255 FORMAT(9X, '-1 +1 -1 +1') c if(itypf.eq.1.and.iwfit.gt.0)then ! freq fit,weighted WRITE(3,258) WRITE(3,1258) WRITE(header1,258) WRITE(header2,1258) endif 258 FORMAT( * 6X,'J K K <- J K K F MF<- F MF',4X,'Volts',8x, * 'Obs',10X,'Obs-Calc',' Error ',' Calc ',6x, * 'Field Field**2') 1258 FORMAT(9X, '-1 +1 -1 +1') c if(itypf.eq.2.and.iwfit.eq.0)then ! diff fit,unweighted WRITE(3,412) WRITE(3,1412) WRITE(header1,412) WRITE(header2,1412) endif 412 FORMAT( * 6X,'J K K <- J K K F MF<-F MF',3X,'Volts Frequency', * 6x,'Stark shift',7x,'Frequency',7x,'Field Field**2') 1412 FORMAT( * 9X, '-1 +1 -1 +1',T50, * 'Obs Obs Obs-Calc Calc') c if(itypf.eq.2.and.iwfit.gt.0)then ! freq fit,weighted WRITE(3,259) WRITE(3,1259) WRITE(header1,259) WRITE(header2,1259) endif 259 FORMAT( * 6X,'J K K <- J K K F MF<-F MF',3X,'Volts Frequency', * 6x,'Stark shift',5x,' Error ',' Frequency',7x, * 'Field Field**2') 1259 FORMAT( * 9X, '-1 +1 -1 +1',T50, * 'Obs Obs Obs-Calc Calc') c write(*,'(1x)') C C...write out frequencies C SUMSQ=0.D0 SUMWSQ=0.d0 nspace=1 ncomm=2 sfmin= 1.D20 sfmax=-1.D20 DO 7 I=1,NT c if(npreds.ge.1.and.i.eq.1)then nsym=1 write(3,'(1x/a,i1/1x)')'Plot symbol: ',nsym if(npreds.ge.2)then maxnc=len_trim(filnam) if(maxnc.gt.27)maxnc=27 do 70 nc=1,maxnc if(filnam(nc:nc).eq.' '.or.filnam(nc:nc).eq.'.')goto 71 70 continue 71 nc=nc-1 if(nc.gt.26)nc=26 filout=filnam(1:nc)//'1'//'.dat' open(2,file=filout,status='unknown') if(npreds.eq.2)write(2,72)coment,(inq(i,j),j=1,10), * 'Volts','Volts' if(npreds.eq.3)write(2,72)coment,(inq(i,j),j=1,10), * 'field','field' 72 format('! ',a/'!',78(1h-)/'!'/ * '! Stark component: ',3i3,' <-',3i3, 2x,2i3,' <-',2i3/ * '!'/'!'/'! Voltage shift calc.freq. ', * a,6x,a,'**2'/'!') write(descdat(nsym),73)(inq(i,j),j=1,10) 73 format( * '! Stark component: ',3i3,' <-',3i3, 2x,2i3,' <-',2i3) endif endif c if(npreds.ge.1.and.i.ne.1.and. * ( (index(i).ne.index(i-1)).or.(indexf(i).ne.indexf(i-1)) * .or.(inq(i,8).ne.inq(i-1,8)) * .or.(inq(i,10).ne.inq(i-1,10))) )then nsym=nsym+1 if(nsym.lt.10)then write(3,'(1x/a,i1/1x)')'Plot symbol: ',nsym if(npreds.ge.2)then write(filout(nc+1:nc+1),'(i1)')nsym filout=filout(1:nc+1)//'.dat' close(2) open(2,file=filout,status='unknown') if(npreds.eq.2)write(2,72)coment,(inq(i,j),j=1,10), * 'Volts','Volts' if(npreds.eq.3)write(2,72)coment,(inq(i,j),j=1,10), * 'field','field' endif else write(3,'(1x/2a/1x)') * 'Plot symbol: ',char(ichar('A')+nsym-10) if(npreds.ge.2)then write(filout(nc+1:nc+2),'(i2)')nsym filout=filout(1:nc+2)//'.dat' close(2) open(2,file=filout,status='unknown') if(npreds.eq.2)write(2,72)coment,(inq(i,j),j=1,10), * 'Volts','Volts' if(npreds.eq.3)write(2,72)coment,(inq(i,j),j=1,10), * 'field','field' endif endif write(descdat(nsym),73)(inq(i,j),j=1,10) endif c if(itypf.ne.2)then TFREQ(I,3)=TFREQ(I,1)-TFREQ(I,2) else TFREQ(I,3)=(TFREQ(I,1)-TFREQ(ifree(i),1))- * (TFREQ(I,2)-TFREQ(ifree(i),2)) error(i)=dsqrt(referr(i)**2+referr(ifree(i))**2) if(error(i).ne.0.d0)then if(i.eq.ifree(i))error(i)=referr(i) wt(i)=1.d0/error(i)**2 else error(i)=0.d0 wt(i)=1.d0 endif endif c if(ispace(nspace).eq.i)then nspace=nspace+1 write(3,'(1x)') endif c 33 if(icomml(ncomm).eq.i)then write(3,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 33 endif c sfield=volts(i)*(1.0+con(36)*volts(i))/con(35) IF(LINFIT(I).NE.1)GOTO 251 c c---fitted lines c SUMSQ=SUMSQ+TFREQ(I,3)**2 SUMWSQ=SUMWSQ+TFREQ(I,3)**2*WT(I) c c...fit to absolute frequencies c if(itypf.eq.1.and.iwfit.eq.0)then ! unweighted WRITE(3,211)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3),tfreq(i,2),sfield,sfield*sfield WRITE(linout,211)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3),tfreq(i,2),sfield,sfield*sfield endif 211 FORMAT(I3,'. ',3I3,' ',3I3,2X,2(I3),' ',2(I3), * F10.1,F15.5,F13.5,f15.4,F12.4,f12.2) c if(itypf.eq.1.and.iwfit.gt.0)then ! weighted WRITE(3,360)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3),error(i),tfreq(i,2),sfield,sfield*sfield WRITE(linout,360)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3),error(i),tfreq(i,2),sfield,sfield*sfield endif 360 FORMAT(I3,'. ',3I3,' ',3I3,2X,2(I3),' ',2(I3), * F10.1,F15.5,F13.5,f12.5,f15.4,F12.4,f12.2) c c...fit to frequency differences c if(itypf.eq.2)then shifto=TFREQ(i,1)-TFREQ(ifree(i),1) c if(iwfit.eq.0)then ! unweighted WRITE(3,414)I,(INQ(I,J),J=1,10),volts(i), * TFREQ(I,1),shifto,TFREQ(I,3),tfreq(i,2), * sfield,sfield*sfield WRITE(linout,414)I,(INQ(I,J),J=1,10),volts(i), * TFREQ(I,1),shifto,TFREQ(I,3),tfreq(i,2), * sfield,sfield*sfield endif 414 FORMAT(I3,'. ',3I3,' ',3I3,2X,4(I3),F8.1,F13.4,2F9.4, * f15.4,f12.4,f12.2) c if(iwfit.gt.0)then ! weighted WRITE(3,361)I,(INQ(I,J),J=1,10),volts(i), * TFREQ(I,1),shifto,TFREQ(I,3),error(i),tfreq(i,2), * sfield,sfield*sfield WRITE(linout,361)I,(INQ(I,J),J=1,10),volts(i), * TFREQ(I,1),shifto,TFREQ(I,3),error(i),tfreq(i,2), * sfield,sfield*sfield endif 361 FORMAT(I3,'. ',3I3,' ',3I3,2X,4(I3),F8.1,F13.4,2F9.4, * f12.4,f15.4,f12.4,f12.2) endif c omc=abs(tfreq(i,3)) do 561 n=1,5 if(omc.gt.worstomc(n))then do 562 nn=6,n+1,-1 worstomc(nn)=worstomc(nn-1) worstlin(nn)=worstlin(nn-1) 562 continue worstomc(n)=omc worstlin(n)=linout goto 566 endif 561 continue c 566 GOTO 7 c c---unfitted lines c 251 if(itypf.eq.1.and.iwfit.eq.0) * WRITE(3,256)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3),tfreq(i,2),sfield,sfield*sfield 256 FORMAT(I3,'. ',3I3,' ',3I3,2X,2(I3),' ',2(I3), * F10.1,F15.5,F13.5,'--',f13.4,f12.4,f12.2) if(itypf.eq.1.and.iwfit.gt.0) * WRITE(3,362)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3),error(i),tfreq(i,2),sfield,sfield*sfield 362 FORMAT(I3,'. ',3I3,' ',3I3,2X,2(I3),' ',2(I3), * F10.1,F15.5,F13.5,'--',f10.5,f15.4,f12.4,f12.2) c if(itypf.eq.1)then if(npreds.eq.2)WRITE(2,2133)volts(i),-TFREQ(I,3), * tfreq(i,2),volts(i),volts(i)*volts(i) if(npreds.eq.3)WRITE(2,2133)volts(i),-TFREQ(I,3), * tfreq(i,2),sfield,sfield*sfield 2133 FORMAT(F10.1,F13.5,f15.4,F12.4,f14.2) endif c if(itypf.eq.2)then shifto=TFREQ(i,1)-TFREQ(ifree(i),1) if(iwfit.eq.0)WRITE(3,415)I,(INQ(I,J),J=1,10),volts(i), * TFREQ(I,1),shifto,TFREQ(I,3),tfreq(i,2), * sfield,sfield*sfield 415 FORMAT(I3,'. ',3I3,' ',3I3,2X,4(I3),F8.1,F13.4,2F9.4,'--', * f13.4,F12.4,f12.2) if(iwfit.gt.0)WRITE(3,363)I,(INQ(I,J),J=1,10),volts(i), * TFREQ(I,1),shifto,TFREQ(I,3),error(i),tfreq(i,2), * sfield,sfield*sfield 363 FORMAT(I3,'. ',3I3,' ',3I3,2X,4(I3),F8.1,F13.4,2F9.4,'--', * f10.4,f15.4,F12.4,f12.2) C if(npreds.eq.2)WRITE(2,2133)volts(i),-TFREQ(I,3), * tfreq(i,2),volts(i),volts(i)*volts(i) if(npreds.eq.3)WRITE(2,2133)volts(i),-TFREQ(I,3), * tfreq(i,2),sfield,sfield*sfield endif c if(tfreq(i,2).gt.sfmax)sfmax=tfreq(i,2) if(tfreq(i,2).lt.sfmin)sfmin=tfreq(i,2) 7 CONTINUE c close(2) c c if(npreds.eq.0)then SDEVSQ=DSQRT(SUMSQ/NDEGF) SDEVW=DSQRT(SUMWSQ/NDEGF) WRITE(3,212)SDEVSQ WRITE(*,212)SDEVSQ if(iwfit.gt.0)then WRITE(3,213)SDEVW WRITE(*,213)SDEVW else write(3,'(1x)') write(*,'(1x)') endif endif 212 FORMAT(1X/' Standard deviation = ',F14.6) 213 FORMAT( ' Weighted deviation = ',F14.6/) C C...Determine increments to constants by linear least squares (works C satisfactorily but could probably be replaced by a nonlinear routine from C Numerical Recipes) C CALL LSQRW(NT,NCON,SD,X) ! <---- LSQRW DO 10 I=1,NCON J=ICFIT(I) if( (j.ne.34) .and. (j.ne.27) .and. (j.ne.28) )then CON(J)=CON(J)+X(I) else CON(J)=CON(J)-X(I) endif 10 CONTINUE C C...Write constants predicted by each iteration C if(npreds.eq.0)then WRITE(3,11)LP WRITE(*,11)LP 11 FORMAT (1X/' ITERATION NO =',I2,' CONSTANTS, deviations and', 1 ' changes:'/) DO 13 I=1,NCON J=ICFIT(I) WRITE(*,12)PO(J),CON(J),SD(I),X(I) WRITE(3,12)PO(J),CON(J),SD(I),X(I) 12 FORMAT (1X,A,' =',F20.9,' +-',F17.9,F22.9) 13 CONTINUE endif c 15 CONTINUE c C--------------------------------------------------------- C C Write final results C C...Constants, with additional conversion of values 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 WRITE(3,16) 16 FORMAT (1X//' FINAL RESULTS OF LEAST SQUARES FITTING PROCEDURE' 1 /1X,48(1H=)//' FITTED CONSTANTS:'/) C L=0 DO 1002 N=1,NCONST UNIT=NAMUN(UNITS(N)) if(units(n).lt.7)then ASCLD=CON(N)*10.D0**(3*(UNITS(N)-1)) else ascld=con(n) endif c IF(IA(N).EQ.0)GOTO 1003 C C...value with error C L=L+1 if(units(n).lt.7)then ESCLD=SD(L)*10.D0**(3*(UNITS(N)-1)) else escld=sd(L) endif WRITE(CONVAL,'(F25.16)')ASCLD WRITE(ERVAL,'(F25.16)')ESCLD c CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR,nerd) ! <---- CONFOR c nchars=17+ndcon+nderor if(nchars.le.39)then WRITE(OUTCON(N),411)POSHRT(N),UNIT,CONVAL(1:NDCON), * ERVAL(1:nderor) else write(outcon(n),'(2a)')conval(1:ndcon),erval(1:nderor) endif GOTO 1002 C C...value without error C 1003 call CONFIX(con(n),ascld,conval,iend) ! <---- CONFIX c WRITE(OUTCON(N),1005)POSHRT(N),UNIT,CONVAL 1002 CONTINUE c c...two column output (blank unused constants) c LLLL=NCONST/2 DO 1107 N=1,LLLL if(outcon(n)(1:6).eq.' ')outcon(n)=' ' if(outcon(n+llll)(1:6).eq.' ')outcon(n+llll)=' ' WRITE(3,1108)OUTCON(N),OUTCON(N+LLLL) 1107 CONTINUE IF(2*LLLL.NE.NCONST)WRITE(3,1108)' ',OUTCON(NCONST) 1108 FORMAT(1X,A38,1X,A38) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C write(3,'(1x)') IF(CON(17).NE.0.D0.OR.CON(18).NE.0.0D0)THEN XCC=-0.5D0*(CON(17)+CON(18)) XBB=CON(18)+XCC WRITE(3,115)'1:',XBB,'1:',XCC ENDIF 115 FORMAT(1X,A,'X.bb /MHz ',F23.12 1 /1X,A,'X.cc /MHz ',F23.12) C C...Correlation coefficients C if(NCON.gt.1)then WRITE(3,121) 121 FORMAT(/' CORRELATION COEFFICIENTS:') NMIN=1 NMAX=8 124 IF(NMAX.GT.NCON)NMAX=NCON WRITE(3,'(1X)') WRITE(3,120)(POSHRT(ICFIT(I)),I=NMIN,NMAX) 120 FORMAT(7X,8(3X,A6)) WRITE(3,'(1X)') DO 122 I=1,NCON LMAX=NMAX IF(I.LT.LMAX)LMAX=I IF(LMAX.LT.NMIN)GOTO 122 WRITE(3,123)POSHRT(ICFIT(I)),(CORELC(I,J),J=NMIN,LMAX) 122 CONTINUE 123 FORMAT(1X,A6,8(F9.4)) IF(NMAX.EQ.NCON)GOTO 125 NMIN=NMIN+8 NMAX=NMAX+8 GOTO 124 endif C C...Worst lines C 125 WRITE(3,564) 564 FORMAT(//' WORST FITTED LINES:'/) write(3,'(a)')header1(1:len_trim(header1)) write(3,'(a)')header2(1:len_trim(header2)) write(3,'(1x)') do 563 n=1,5 write(3,'(a)')worstlin(n)(1:len_trim(worstlin(n))) 563 continue C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C...Diagnostic plots for the predictive mode C if(npreds.ge.1)then c volmin=50000. volmax=0. fmax=-1.D20 fmin= 1.D20 fzmin= 1.d20 fzmax=-1.d20 do 50 i=1,nt if(volts(i).lt.volmin)volmin=volts(i) if(volts(i).gt.volmax)volmax=volts(i) if(-tfreq(i,3).lt.fmin)fmin=-tfreq(i,3) if(-tfreq(i,3).gt.fmax)fmax=-tfreq(i,3) if(volts(i).lt.1.0d0)then if(tfreq(i,2).lt.fzmin)fzmin=tfreq(i,2) if(tfreq(i,2).gt.fzmax)fzmax=tfreq(i,2) endif 50 continue if(flower.ne.0.d0.or.fupper.ne.0.0d0)then if(flower.lt.fupper)then fmin=flower fmax=fupper endif endif c c...header c plotbf(1)=' ' plotbf(nhpix)=' ' do 57 i=2,nhpix-1 plotbf(i)='-' 57 continue write(3,'(1x//132a1)')(plotbf(j),j=1,nhpix) do 68 i=2,nhpix-1 plotbf(i)=' ' 68 continue write(comcop,'(''| '',a)')coment plotbf(nhpix)='|' write(3,'(132a1)')(plotbf(j),j=1,nhpix) do 67 i=2,nhpix-1 plotbf(i)=' ' 67 continue write(3,'(132a1)')(plotbf(j),j=1,nhpix) write(cbuff,56) '| Voltage (X): ',volmin,volmax do 58 i=51,nhpix-1 plotbf(i)=' ' 58 continue plotbf(nhpix)='|' if(itypre.eq.1)then write(clabel,'(a)')'LINEAR SCALE ' else write(clabel,'(a)')'QUADRATIC SCALE' endif write(3,'(132a1)')(plotbf(j),j=1,nhpix) write(clabel,'(a)')' ' write(cbuff,56) '| Field : ', * volmin/con(35),volmax/con(35) write(3,'(132a1)')(plotbf(j),j=1,nhpix) write(cbuff,56) '| Stark shift (Y): ', fmin,fmax write(3,'(132a1)')(plotbf(j),j=1,nhpix) 56 format(a,2f15.4) c c...body of plot c fincr=(fmax-fmin)/real(nvpix-1) flow=fmax-0.5d0*fincr fhigh=fmax+0.5d0*fincr c c...determine vertical markers vmark=0.001d0 nmult=0 60 nmark=int( (volmax-volmin)/vmark ) if(nmark.gt.5)then vmark=vmark*2.d0 nmult=nmult+1 if(nmult.eq.2)then vmark=vmark*1.25d0 nmult=-1 endif goto 60 endif write(cbuff,61) '| Voltage marker: ', vmark write(3,'(132a1)')(plotbf(j),j=1,nhpix) c c...determine horizontal markers fmark=0.001d0 nmult=0 63 nmark=int( (fmax-fmin)/fmark ) if(nmark.gt.5)then fmark=fmark*2.d0 nmult=nmult+1 if(nmult.eq.2)then fmark=fmark*1.25d0 nmult=-1 endif goto 63 endif write(cbuff,61) '| Frequency marker: ', fmark write(3,'(132a1)')(plotbf(j),j=1,nhpix) 61 format(a,15x,f15.4) fcurm=int(fmax/fmark)*fmark if(fcurm.ge.fmax)fcurm=fcurm-fmark c c...box and grid do 51 i=1,nvpix do 53 j=1,nhpix plotbf(j)=' ' 53 continue if(i.eq.1.or.i.eq.nvpix)then do 55 j=1,nhpix plotbf(j)='-' 55 continue endif if(fcurm.le.fhigh.and.fcurm.gt.flow)then if(i.ne.1.and.i.ne.nvpix)then do 62 j=2,nhpix-1 plotbf(j)='.' 62 continue fcurm=fcurm-fmark endif endif if(0.d0.le.fhigh.and.0.d0.gt.flow)then do 54 j=2,nhpix-1 plotbf(j)=',' 54 continue endif c if(i.eq.1.or.i.eq.nvpix)goto 65 vcurm=int(volmin/vmark)*vmark if(vcurm.le.volmin)vcurm=vcurm+vmark 64 if(vcurm.ge.volmax)goto 65 if(itypre.eq.1)then ibuf=int( (vcurm-volmin)*nhpix/(volmax-volmin) ) else ibuf=int((vcurm**2-volmin**2)*nhpix/(volmax**2-volmin**2)) endif plotbf(ibuf)='.' vcurm=vcurm+vmark goto 64 65 plotbf(1)='|' plotbf(nhpix)='|' c nsym=1 do 52 j=1,nt c c...increment data set marker if(npreds.ge.1.and.j.ne.1.and. * ( (index(j).ne.index(j-1)).or.(indexf(j).ne.indexf(j-1)) * .or.(inq(j,8).ne.inq(j-1,8)) * .or.(inq(j,10).ne.inq(j-1,10))) )nsym=nsym+1 c c...fixed scaling: mark +ve overflow on the top margin if(flower.ne.fupper.and.i.eq.1)then if(-tfreq(j,3).gt.fhigh)then if(itypre.eq.1)then ibuf=int( (volts(j)-volmin)*nhpix/(volmax-volmin) ) if(ibuf.lt.1)ibuf=1 if(ibuf.gt.nhpix)ibuf=nhpix else ibuf=int( (volts(j)**2-volmin**2)*nhpix/ * (volmax**2-volmin**2) ) if(ibuf.lt.1)ibuf=1 if(ibuf.gt.nhpix)ibuf=nhpix endif if(nsym.lt.10)then write(plotbf(ibuf),'(i1)')nsym else write(plotbf(ibuf),'(a1)')char(ichar('A')+nsym-10) endif endif endif c c...fixed scaling: mark -ve overflow on the bottom margin if(flower.ne.fupper.and.i.eq.nvpix)then if(-tfreq(j,3).le.flow)then if(itypre.eq.1)then ibuf=int( (volts(j)-volmin)*nhpix/(volmax-volmin) ) if(ibuf.lt.1)ibuf=1 if(ibuf.gt.nhpix)ibuf=nhpix else ibuf=int( (volts(j)**2-volmin**2)*nhpix/ * (volmax**2-volmin**2) ) if(ibuf.lt.1)ibuf=1 if(ibuf.gt.nhpix)ibuf=nhpix endif if(nsym.lt.10)then write(plotbf(ibuf),'(i1)')nsym else write(plotbf(ibuf),'(a1)')char(ichar('A')+nsym-10) endif endif endif c c...mark data points in the plot box if(-tfreq(j,3).le.fhigh.and.-tfreq(j,3).gt.flow)then if(itypre.eq.1)then ibuf=int( (volts(j)-volmin)*nhpix/(volmax-volmin) ) if(ibuf.lt.1)ibuf=1 if(ibuf.gt.nhpix)ibuf=nhpix else ibuf=int( (volts(j)**2-volmin**2)*nhpix/ * (volmax**2-volmin**2) ) if(ibuf.lt.1)ibuf=1 if(ibuf.gt.nhpix)ibuf=nhpix endif if(nsym.lt.10)then write(plotbf(ibuf),'(i1)')nsym else write(plotbf(ibuf),'(a1)')char(ichar('A')+nsym-10) endif endif 52 continue write(3,'(132a1)')(plotbf(j),j=1,nhpix) flow=flow-fincr fhigh=fhigh-fincr 51 continue c c...gle output (NPREDS=2,3) if(npreds.ge.2)then eldist=con(35) if(npreds.eq.2)eldist=1.d0 fzmean=(fzmin+fzmax)*0.5d0 write(*,290)fzmin,fzmax,fzmean 290 format(' Zero field frequencies:',T30,'minimum =',F12.4/ * T30,'maximum =',F12.4/ * T30,'mean =',F12.4) call gleout(filout,nc,nsym,volmin,volmax,eldist,sfmin,sfmax, ! <---- gleout * fmin,fmax,fzmean,flower,fupper,vmark,fmark,itypre,iveho) endif c c endif C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C C...Write the EXP.DAT file (if in fitting mode) C if(npreds.eq.0)then OPEN(10,FILE='EXP.DAT',STATUS='UNKNOWN') write(10,1123)filnam(1:LEN_TRIM(FILNAM)) 1123 format('! Experimental points from file: ',a/'!'/ * '! voltage frequency'/'!') c do 1122 i=1,nt if(linfit(i).eq.0)goto 1122 write(10,'(f12.2,f16.4)')volts(i),tfreq(i,1) 1122 continue c close(10) endif C C...Write also the EXPTPLOT.DAT file (if in fitting mode) C if(npreds.eq.0)then OPEN(10,FILE='EXPTPLOT.DAT',STATUS='UNKNOWN') write(10,30)filnam(1:LEN_TRIM(FILNAM)) 30 format('! Experimental points from file: ',a/ * '!'/'! voltage deltaF F ', * 'E E**2'/'!') c eldist=con(35) c if(npreds.eq.2)eldist=1.d0 eldist=1.0d0 do 2122 i=1,nt if(linfit(i).eq.0)goto 2122 v=volts(i) obsf=tfreq(i,1) if(v.eq.0.d0)fzero=obsf if(eldist.eq.1.d0)then write(10,2135)v,obsf-fzero,obsf,v,v*v else sfield=v/eldist write(10,2135)v,obsf-fzero,obsf,sfield,sfield*sfield endif 2135 FORMAT(F10.1,F13.5,f15.4,F12.4,f14.2) 2122 continue c close(10) endif C C C...Contributions to frequencies (optional) C IF(ICONTR.LT.1)GOTO 126 C DO 1130 J=1,NCON ! Fortran2018 K=ICFIT(J) DO 130 I=1,NT DER(I,J)=DER(I,J)*CON(K) if(k.ge.32.and.k.le.34)der(i,j)=der(i,j)*0.5d0 if(k.eq.35)der(i,j)=-der(i,j)*0.5d0 if(k.eq.22.or.k.eq.19)der(i,j)=der(i,j)*0.5d0 130 CONTINUE 1130 CONTINUE ! Fortran2018 C IF(ICONTR.EQ.2)THEN WRITE(3,153) 153 FORMAT(1X/' TRANSITIONS:'// 1 7X,'J K K <- J K K I F<- I F',4X,'Volts',8x, 2 'Obs',8X,'Field free'/ 3 10X, '-1 +1 -1 +1'/) DO 151 I=1,NT TFREQ(I,3)=TFREQ(I,1) DO 152 J=1,NCON K=ICFIT(J) IF(k.GE.32)TFREQ(I,3)=TFREQ(I,3)-DER(I,J) 152 CONTINUE WRITE(3,311)I,(INQ(I,J),J=1,10),volts(i),TFREQ(I,1), * TFREQ(I,3) 151 CONTINUE GOTO 126 ENDIF 311 FORMAT(I3,'. ',3I3,' ',3I3,2X,2(I3),' ',2(I3), 1 F10.1,2F14.4) C WRITE(3,131) 131 FORMAT(//' CONTRIBUTIONS (from gradients) OF INDIVIDUAL ', * 'CONSTANTS TO FREQUENCY:'/) NMIN=1 NMAX=8 132 IF(NMAX.GT.NCON)NMAX=NCON WRITE(3,150)(POSHRT(ICFIT(I)),I=NMIN,NMAX) 150 FORMAT(67X,8(3X,A6)) IF(NMAX.EQ.NCON)GOTO 133 NMIN=NMIN+8 NMAX=NMAX+8 GOTO 132 133 WRITE(3,'(1X)') C nspace=1 ncomm=2 DO 134 N=1,NT FORMAL=F1 NMIN=1 NMAX=8 138 IF(NMAX.GT.NCON)NMAX=NCON DO 135 L=NMIN,NMAX LLL=5 DERIV=DABS(DER(N,L)) DIV=10.D0 137 IF(DERIV/DIV.LT.1.D0)GOTO 136 DIV=DIV*10.D0 LLL=LLL-1 GOTO 137 136 IF(LLL.LT.0)LLL=0 LLLL=L IF(LLLL.GT.8)LLLL=LLLL-8 FORMAL(55+LLLL*5:55+LLLL*5)=CHAR(LLL+ICHAR('0')) 135 CONTINUE IF(NMIN.GT.1)GOTO 139 c if(ispace(nspace).eq.n)then nspace=nspace+1 write(3,'(1x)') endif c 313 if(icomml(ncomm).eq.n)then write(3,'(a)')comblk(icomm(ncomm-1)+1:icomm(ncomm)) ncomm=ncomm+1 goto 313 endif c WRITE(linout,FORMAL)N,(INQ(N,J),J=1,10),volts(N),TFREQ(N,1), * (DER(N,L),L=NMIN,NMAX) if(linfit(n).eq.0)linout(64:65)='--' WRITE(3,'(a)')linout(1:lentrm(linout)) c GOTO 140 139 WRITE(3,FORMAL)(DER(N,L),L=NMIN,NMAX) C 140 IF(NMAX.EQ.NCON)GOTO 134 FORMAL=F2 NMIN=NMIN+8 NMAX=NMAX+8 GOTO 138 134 CONTINUE C C 126 if(npreds.lt.1)write(3,'(1x/1x,78(1h-))') CLOSE(3) C Write(*,'(1X//)') C STOP END C C_____________________________________________________________________________ C SUBROUTINE gleout(filout,nc,nsym,volmin,volmax,eldist,sfmin,sfmax, * fmin,fmax,fzmean,flower,fupper,vmark,fmark,itypre,iveho) C C filout - on input carries the name of the last .DAT file c nc - length of generic file name c nsym - number of Stark components C volmin,volmax - voltage limits c eldist - electrode separation (if set to 1 then plot of c frequency versus voltage) c sfmin,sfmax - absolute frequency limits of points in data c fmin,fmax - frequency shift limits c fzmean - mean zero field frequency for calculated Stark components c flower,fupper - predefined frequency shift limits c vmark,fmark - voltage and frequency markers c itypre - type of plot =1 linear =2 quadratic c iveho - horizontal (=0) or vertical (=1) gle plot C implicit real(8) (a-h,o-z) character filout*30,cns*2,line*120,linea*120,lineb*120, * linev*120,coment*78 chARACTER descdat*59(200) c common /cmnt/coment,descdat c c...change horizontal limits and marker for field plot c if(eldist.ne.1.d0)then volmin=volmin/eldist volmax=volmax/eldist vmark=vmark/eldist vmark=0.001d0 nmult=0 60 nmark=int( (volmax-volmin)/vmark ) if(nmark.gt.5)then vmark=vmark*2.d0 nmult=nmult+1 if(nmult.eq.2)then vmark=vmark*1.25d0 nmult=-1 endif goto 60 endif endif c 101 write(*,100,advance='NO') 100 format(1x/' Specify type of frequency axis:'// * 10x,'0 = difference from zero field frequency (Stark shift)'/ * 10x,'1 = absolute frequency'// * 6x,'... ') read(*,'(i5)',err=101)iftype if(iftype.lt.0.or.iftype.gt.1)goto 101 c filout=filout(1:nc)//'.gle' open(2,file=filout,status='unknown') write(2,102)coment 102 format('!',79(1H-)/'! Script generated by QSTARK for:'/'!' * /'! ',a/'!',79(1H-)/'!') c if(iveho.eq.0)then write(2,2)'size 29.5 21' else write(2,2)'size 21 29.5' endif 2 format(a) 3 format(a/) 4 format(4x,a) write(2,2)'lw=0.03' write(2,2)'sizem=0.4' write(2,2)'tlen=0.5' write(2,2)'stlen=tlen*0.5' write(2,3)'set lwidth lw' write(2,3)'set hei 0.6' if(iveho.eq.0)then write(2,2)'amove 0.5 0.5' else write(2,2)'amove 1.0 1.0' endif write(2,2)'begin graph' c write(2,4)'nobox' if(iveho.eq.0)then write(2,4)'size 31 22' else write(2,4)'size 22 31' endif do 1 i=1,nsym write(cns,'(i2)')i if(i.lt.10)then filout=filout(1:nc)//cns(2:2)//'.dat' if(itypre.eq.1)then if(iftype.eq.0)write(2,4) * 'data '//filout(1:nc+5)//' d'//cns(2:2)//'=c4,c2' if(iftype.eq.1)write(2,4) * 'data '//filout(1:nc+5)//' d'//cns(2:2)//'=c4,c3' else if(iftype.eq.0)write(2,4) * 'data '//filout(1:nc+5)//' d'//cns(2:2)//'=c5,c2' if(iftype.eq.1)write(2,4) * 'data '//filout(1:nc+5)//' d'//cns(2:2)//'=c5,c3' endif else filout=filout(1:nc)//cns(1:2)//'.dat' if(itypre.eq.1)then if(iftype.eq.0)write(2,4) * 'data '//filout(1:nc+6)//' d'//cns(1:2)//'=c4,c2' if(iftype.eq.1)write(2,4) * 'data '//filout(1:nc+6)//' d'//cns(1:2)//'=c4,c3' else if(iftype.eq.0)write(2,4) * 'data '//filout(1:nc+6)//' d'//cns(1:2)//'=c5,c2' if(iftype.eq.1)write(2,4) * 'data '//filout(1:nc+6)//' d'//cns(1:2)//'=c5,c3' endif endif 1 continue c write(2,2)' ' write(2,4)'ticks color grey40' write(2,4)'dsubticks on' if(itypre.eq.1)then write(line,7)'xaxis min',volmin,' max',volmax,' dticks', * vmark,' dsubticks',vmark*0.1d0 write(2,4)line(1:86)//' grid' if(eldist.eq.1.d0)then write(2,4)'xtitle "\rm Voltage\,/V" '// * 'hei 0.7 dist 0.4 font texcmmi' else write(2,4)'xtitle "E\,\rm/Vcm^{-1}" '// * 'hei 0.7 dist 0.4 font texcmmi' endif write(2,4)'xticks length tlen' write(2,4)'xsubticks length stlen' else write(line,7)'xaxis min',volmin**2,' max',volmax**2 * ,' dticks', * 4.*volmax**2,' dsubticks',2.*volmax**2 write(2,4)line(1:86)//' grid' if(eldist.eq.1.d0)then write(2,4)'xtitle "\rm Voltage^2\,/V^2" '// * 'hei 0.7 dist 0.4 font texcmmi' else write(2,4)'xtitle "E\rm^2\,/V^2cm^{-2}" '// * 'hei 0.7 dist 0.4 font texcmmi' endif write(2,4)'xticks length tlen' write(2,4)'xlabels dist 0.4' c volm=volmin+vmark line ='xplaces ' linea='xnames ' lc=8 11 if(volm.lt.volmax)then write(line(lc+1:lc+11),'(i11)')nint(volm**2) write(linea(lc+1:lc+11),'(a,i5,a)')' "',nint(volm),'^2 "' write(lineb(lc+1:lc+11),'(a,i5,a)')' " " ' volm=volm+vmark lc=lc+11 goto 11 endif write(2,4)line(1:lc) write(2,4)linea(1:lc) line(1:8) ='x2places' lineb(1:8)='x2names ' write(2,4)line(1:lc) write(2,4)lineb(1:lc) endif 7 format(a,f14.1,a,f14.1,a,f14.1,a,f14.1) c write(2,2)' ' if(iftype.eq.0)then write(line,6)'yaxis min',fmin,' max',fmax,' dticks',fmark, * ' dsubticks',fmark*0.1d0 else if((flower.ne.0.d0.or.fupper.ne.0.0d0) * .and. (flower.lt.fupper) )then fmin=fzmean+flower fmax=fzmean+fupper if(fmin.gt.fmax)then fmin=fzmean+fupper fmax=fzmean+flower endif else fmin=sfmin-0.05*(sfmax-sfmin) fmax=sfmax+0.05*(sfmax-sfmin) endif c if((fmax-fmin)/fmark.gt.10.d0)then write(line,6)'yaxis min',fmin,' max',fmax else write(line,6)'yaxis min',fmin,' max',fmax,' dticks',fmark, * ' dsubticks',fmark*0.1d0 endif endif 6 format(a,f10.3,a,f10.3,a,f10.3,a,f10.3) write(2,4)line(1:70)//' grid' if(iftype.eq.0)then write(2,4)'ytitle "\Delta\nu_s\, \rm /MHz" '// * 'hei 0.7 dist 0.4 font texcmmi' else write(2,4)'ytitle "\nu_s\, \rm /MHz" '// * 'hei 0.7 dist 0.4 font texcmmi' endif write(2,4)'yticks length tlen' write(2,4)'ysubticks length stlen' c write(2,2)'!' do 10 i=1,nsym write(cns,'(i2)')i if (i.lt.10)then write(2,4)'d'//'0'//cns(2:2)//' lstyle 1 lwidth lw smooth'// * ' color black '//descdat(i) else write(2,4)'d'//cns(1:2)//' lstyle 1 lwidth lw smooth'// * ' color black '//descdat(i) endif 10 continue c c...produce the EXPTPLOT.DAT file from EXP.DAT (if available) C line(1:5)='! ' filout='xxxxx.xx' open(7,file='exp.dat',status='old',err=21) c npts=0 fshift=0.d0 c c__top of loop reading and processing EXP.DAT c 20 read(7,'(a)',err=1020,end=22)linev c if(linev(1:1).eq.'!')goto 20 if(len_trim(linev).lt.10)goto 20 c read(linev,*,err=1020,end=22)v,f c if(v.eq.0.)fzero=f npts=npts+1 if(npts.eq.1)then open(1,file='exptplot.dat',status='unknown',err=22) write(1,30) 30 format('!'/'! voltage deltaF F ', * 'E E**2'/'!') line(1:1)=' ' endif if(eldist.eq.1.d0)then write(1,2133)v,f-fzero,f,v,v*v else sfield=v/eldist write(1,2133)v,f-fzero,f,sfield,sfield*sfield endif goto 20 2133 FORMAT(F10.1,F13.5,f15.4,F12.4,f14.2) c 1020 write(*, * '(1x//'' WARNING: Input error from EXP.DAT at line: ''/1x,a/)') * line(1:len_trim(line)) C 22 close(7) close(1) write(filout,'(f15.5)')fshift 23 if(filout(1:1).eq.' ')then filout=filout(2:) goto 23 endif c c...batch of commented lines which place points on the plot or can be c adapted to do so c 21 write(2,2)'!' write(2,2)'! Lines below allow plotting of experimental points'// * ' on top of' write(2,2)'! calculated lines, uncomment/modify as necessary' write(2,2)'!' if(itypre.eq.1)then if(iftype.eq.0)write(2,2) * line(1:5)//'data '//'exptplot.dat'//' d98=c4,c2' if(iftype.eq.1)write(2,2) * line(1:5)//'data '//'exptplot.dat'//' d98=c4,c3' else if(iftype.eq.0)write(2,2) * line(1:5)//'data '//'exptplot.dat'//' d98=c5,c2' if(iftype.eq.1)write(2,2) * line(1:5)//'data '//'exptplot.dat'//' d98=c5,c3' endif write(2,2)line(1:5)//'let d99 = d98-'//filout write(2,2)line(1:5)//'d99 marker wcircle msize sizem' c write(2,2)'!' write(2,2)'! d1 marker wcircle msize sizem' write(2,2)'! d1 marker fcircle msize sizem' write(2,2)'! d1 marker wsquare msize sizem' write(2,2)'! d1 marker fsquare msize sizem' write(2,2)'! d1 marker wtriangle msize sizem' write(2,2)'! d1 marker ftriangle msize sizem' write(2,2)' ' write(2,2)'end graph' c c...second plot for short ticks for quadratic field c if(itypre.eq.2)then write(2,2)' ' write(2,2)'!' write(2,2)'! Second, overlaid graph for short ticks for'// * ' quadratic field/voltage ' write(2,2)'!' write(2,2)'begin graph' write(2,4)'nobox' if(iveho.eq.0)then write(2,4)'size 31 22' else write(2,4)'size 22 31' endif write(line,7)'xaxis min',volmin**2,' max',volmax**2 * ,' dticks', * 4.*volmax**2,' dsubticks',2.*volmax**2 write(2,4)line(1:86) write(2,4)'xticks length stlen' write(2,4)'xlabels off' c vmark=vmark/5.d0 volm=volmin+vmark 14 line ='xplaces ' linea='xnames ' lc=8 12 if(volm.lt.volmax)then write(line(lc+1:lc+11),'(i11)')nint(volm**2) write(linea(lc+1:lc+11),'(a,i5,a)')' " " ' write(lineb(lc+1:lc+11),'(a,i5,a)')' " " ' volm=volm+vmark lc=lc+11 if(lc.gt.80)goto 15 goto 12 endif 15 write(2,4)line(1:lc) write(2,4)linea(1:lc) c line(1:8) ='x2places' lineb(1:8)='x2names ' write(2,4)line(1:lc) write(2,4)lineb(1:lc) if(volm.lt.volmax)goto 14 c write(2,4)'yaxis off' write(2,2)'end graph' endif c write(2,2)'!' write(2,2)'amove 15. 27.8' write(2,2)'set hei 0.7' write(2,2)'set font texcmr' write(2,2)'text'//descdat(1)(23:30)//' \leftarrow\, '// * descdat(1)(35:42) c write(2,103) 103 format('!'/'!',79(1h-)/'!',79(1h-)) close(2) c RETURN END C_____________________________________________________________________________ C SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR,NERD) C C Constant and error formating for output c c CONVAL - string containing the constant value, this is to be c input in F format and will be replaced on output by result c string of length extending to the last digit of the error c ERVAL - string containing the error value, this is to be c input in F format and will be replaced on output by result c string, this does not contain the decimal point and is meant c to be included in brackets c NDCON - number of digits in the constant value c NDEROR - number of digits in the error value c both NDCON and NDEROR are generated on output c NERD - number of error digits which is to be set on input C c IMPLICIT INTEGER(2) (I-N) CHARACTER(25) CONVAL,ERVAL,CONOUT,EROUT C 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 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) ERVAL(1:NDEROR)=EROUT(1:NDEROR) C RETURN END C C_____________________________________________________________________________ c subroutine confix(const,ascld,confor,iend) c C Fixed constant formatting for output C C This routine is adapted from the PIFORM program such that C 1/ the leading and following square brackets are disabled and C 2/ internal string processing is made in CONVAL but output is to C the shorter string CONFOR as already used by QS c c CONST - unscaled value of the constant c ASCLD - value of the constant scaled by customary multiplier c CONFOR - string for outputting the formatted constant c IEND - length of filled CONVAL, inclusive of the concluding square bracket implicit real(8) (a-h,o-z) character conval*27,confor*25 C if(abs(const).le.1.d-17)ascld=0.d0 IF(ASCLD.NE.0.D0)ASCLD=ASCLD+MOD(ASCLD,1.D-15*ASCLD) WRITE(CONVAL,'(F27.16)')ASCLD c c...deal with leading zero IF(ASCLD.EQ.0.0D0)CONVAL=' 0.0 ' IF(CONVAL(10:10).EQ.' ')CONVAL(10:10)='0' IF(CONVAL(10:10).EQ.'-')THEN CONVAL(10:10)='0' CONVAL(9:9)='-' ENDIF c c...set leading square bracket do 2115 m=1,27 if(CONVAL(M:M).NE.' ')GOTO 2116 2115 continue 2116 if(m.gt.8)then c conval(8:8)='[' conval(8:8)=' ' else if(m.gt.1.and.m.le.8)then c conval(m-1:m-1)='[' conval(m-1:m-1)=' ' endif c c...blank at the end so that number processed is within 15 digit accuracy if(ascld.ne.0.d0)then mmm=int( alog10(real(abs(ascld))) ) if(mmm.ge.-1)mmm=mmm+1 do 2117 m=27,27-mmm,-1 conval(m:m)='0' 2117 continue endif c c...find last non-zero digit DO 1114 M=27,1,-1 IF(CONVAL(M:M).NE.'0'.AND.CONVAL(M:M).NE.' ')GOTO 1115 IF(CONVAL(M:M).EQ.'0')CONVAL(M:M)=' ' 1114 CONTINUE C C...terminating square bracket C 1115 iend=m+1 if(iend.gt.27)iend=27 c conval(iend:iend)=']' conval(iend:iend)=' ' c confor=conval(3:27) c return end C_____________________________________________________________________________ C SUBROUTINE READ(CON,NCON,I1,I2,IDUMP,nt,nqover) C C Routine to: C C - set up the Hamiltonian for given quantum numbers I1,I2 C and spectroscopic constants in CON (the constants that will C actually be used are transferred to CONST) C - find the energy levels required for calculated frequencies C - evaluate the derivatives of the energy levels, and hence the C transitions with respect to the NCON constants to be fitted C c CON - current constants c NCON - number of fitted constants c IDUMP - dump options C NT - number of transitions C C NELINF(i,j) - information on eigenvalues of H required by the spectral data: C (i,1) is the number of transition the level belongs to (-ve for C the lower level) C (i,2) is the number of this eigenvalue among the eigenvalues of H C NEXTR - number of elements of NELINF that are set for each C diagonalization C C NCOL - the number of columns of H that the printer being used can C print on one line C 17 is suitable for a 15" printer in condensed mode C 12 is suitable for 10" in condensed elite (160 characters) c c I1,I2,FH - all in units of 1/2 c C use globdat C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=36,NLIN=2000,NLIN2=2*NLIN,msort=10000) COMMON /BLKEN/TFREQ,DER,LABO,INFOE,volinf,ICFIT,NELEV COMMON /SORTCC/WK(Msort),IPT(msort) common /stark/volts(nlin) DIMENSION LABO(NLIN,2,5), * ICFIT(NCONST),CON(NCONST),CONST(NCONST),TFREQ(NLIN,3), * INFOE(NLIN2,2),DER(NLIN,NCONST),volinf(nlin2) c NCOLS=40 C C C...Set up and diagonalize the Hamiltonian by iterating over required values C of the electric field and, for each, over required values of M.F: in C ascending order from the lowest in the data set to the highest. C The necessary eigenvalues and their derivatives are extracted C following each diagonalization. C istart=1 c c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c MAIN LOOP for each MFCALC and VCALC c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c c...identify the batch of energy levels to be calculated, the lines c commented out allow setting up H with several M.F blocks if necessary c 502 lab=istart mfcalc=infoe(lab,2) C if(con(33).ne.0.d0.or.con(34).ne.0.d0)then C mfmin=mfcalc-2 C if(mfmin.lt.i1-(i1/2)*2)mfmin=i1-(i1/2)*2 C mfmax=mfcalc+2 C else mfmin=mfcalc mfmax=mfcalc C endif vcalc=volinf(lab) jmax=0 c 501 if(infoe(lab,2).ne.mfcalc)goto 500 if(volinf(lab).ne.vcalc)goto 500 labul=2 if(infoe(lab,1).lt.0)labul=1 labpt=iabs(infoe(lab,1)) if(jmax.lt.labo(labpt,labul,2))jmax=labo(labpt,labul,2) lab=lab+1 if(lab.le.nelev)goto 501 c 500 ifin=lab-1 c c...current calculation for VCALC, MFCALC will deliver energy levels from c ISTART to IFIN in INFOE c c...determine the required F limits for the Hamiltonian - calculations will be c made from the lowest value for given M.F to the value resulting from use c of NQOVER c flow=10000 fhigh=mfcalc jlow=10000 jhigh=0 do 1447 ll=istart,ifin NLEV=2 IF(INFOE(LL,1).LT.0)NLEV=1 NLINE=IABS(INFOE(LL,1)) J=LABO(NLINE,NLEV,2) F=labo(nline,nlev,4) if(f.gt.fhigh)fhigh=f if(f.lt.flow)flow=f if(j.gt.jhigh)jhigh=j if(j.lt.jlow)jlow=j 1447 continue fhigh=fhigh+nqover*2 ffhigh=2*jhigh+2*nqover+i1 fhigh=max(fhigh,ffhigh) c flow=flow-nqover*2 1449 if(flow.lt.0)then flow=flow+2 goto 1449 endif fflow=2*jlow-2*nqover-i1 1450 if(fflow.lt.0)then fflow=fflow+2 goto 1450 endif flow=min(flow,fflow) flow=max(flow,mfcalc) c if(idump.ne.0)write(4,1445)mfcalc,vcalc 1445 format(1x/1x,79('-')/' M.F =',i3,'/2',f10.1,' V :') c mstot=0 do 1130 mfloop=mfmin,mfmax,2 ! Fortran2018 do 130 F=flow,fhigh,2 ms=mhigh(i1,i2,f) mstot=mstot+ms if(idump.ge.1)then WRITE(*,221)mfloop,F,ms,ms,MStot,MStot WRITE(4,221)mfloop,F,ms,ms,MStot,MStot endif 130 continue 1130 continue ! Fortran2018 221 FORMAT(10X,'M.F=',I2,'/2',5x,'F=',I2,'/2',5X, * 'block size = ',I3,'x',I3,5x,'H size = ',I4,'x',I4) C...set up and diagonalize H for given values of field and M.F C (decide whether the H matrix will be real, imag=0, or imaginary, imag=1) C DO 105 I=1,NCONST CONST(I)=CON(I) 105 CONTINUE if( (con(34).ne.0.d0) .or. ! mu_c * (con(27).ne.0.d0) .or. ! chi_bc * (con(28).ne.0.d0)) then ! chi_ac imag=1 if(idump.ne.0)then write(*,222)2*mstot,2*mstot ! dump write(4,222)2*mstot,2*mstot ! dump 222 format(48x,' H_augmented size = ',I4,'x',I4) ! dump endif else imag=0 endif c c...Allocate appropriately sized arrays if necessary c IF( 2*MStot.NE.MSIZE )then c if(msize.gt.0)DEALLOCATE(H1,T,NELINF,stat=ialloc) MDIM=MSTOT MSIZE=2*mstot ALLOCATE(H1(MSIZE+5,MSIZE+5),T(MSIZE+5,MSIZE+5), * NELINF(MDIM,2),stat=ialloc) endif c CALL SETHAM(MStot,I1,FLOW,FHIGH,CONST,vcalc,mfmin,mfmax) ! <---- SETHAM c if(idump.ne.0)then write(4,'(1x)') ! dump write(4,1444)(mfloop,mfloop=mfmin,mfmax,2) ! dump write(4,1448)vcalc ! dump 1444 format(10x,'2M.F = ',10I6) ! dump 1448 Format(10x,' V = ',f11.4) ! dump write(4,'(1x)') ! dump if(idump.lt.3)goto 1014 ! dump c c...dump real block c nmin=1 ! dump nmax=mstot ! dump if(imag.ne.0)write(4,'('' A block of H_augmented'')') ! dump if(nmax.gt.ncols)nmax=ncols ! dump 1011 do 1001 i=1,mstot ! dump write(4,1002)i,(h1(i,j),j=nmin,nmax) ! dump 1001 continue ! dump 1002 format( i4,':',60f9.1) ! dump if(nmax.eq.mstot)goto 1010 ! dump nmin=nmin+ncols ! dump nmax=nmax+ncols ! dump if(nmax.gt.mstot)nmax=mstot ! dump write(4,'(1x)') ! dump goto 1011 ! dump 1010 if(imag.ne.1)goto 1014 ! dump c c...dump imaginary block (if there is one) c write(4,'('' B block of H_augmented'')') ! dump nmin=mstot+1 ! dump nmax=mstot+mstot ! dump if(nmax-mstot.gt.ncols)nmax=mstot+ncols ! dump 1072 do 1071 i=1,mstot ! dump write(4,1002)i,(h1(i,j),j=nmin,nmax) ! dump 1071 continue ! dump if(nmax.eq.mstot+mstot)goto 1070 ! dump nmin=nmin+ncols ! dump nmax=nmax+ncols ! dump if(nmax.gt.mstot+mstot)nmax=mstot+mstot ! dump write(4,'(1x)') ! dump goto 1072 ! dump 1070 continue ! dump endif c 1014 IEGEN=0 IORD=0 if(imag.eq.1)then ndiag=2*mstot else ndiag=mstot endif CALL HDIAGM(ndiag,IEGEN,IORD) ! <---- HDIAGM write(*,'('' finished '',i4,''x'',i4,'' M='',i2, * ''/2'',F10.1,'' V''/)')ndiag,ndiag,mfcalc,vcalc c if(idump.ne.0)then ! dump write(4,'('' RAW E-levels:'')') ! dump write(4,'(i4,'':'',6f12.3)') ! dump * (nn,(h1(nnn,nnn),nnn=nn,nn+5),nn=1,mstot,6) ! dump write(4,'(1x)') ! dump if(imag.eq.1)then write(4,'('' RAW E-levels (second set from H_augm):'')') ! dump write(4,'(i4,'':'',6f12.3)') ! dump * (nn,(h1(nnn,nnn),nnn=nn,nn+5),nn=mstot+1,mstot+mstot,6) ! dump write(4,'(1x)') ! dump endif endif ! dump c if(idump.ge.2)then ! dump write(4,'('' Pointers from eigenvectors:'')') ! dump do 1031 nn=1,mstot nnn=1 do 1040 nnnn=2,mstot if(t(nn,nnnn)**2.gt.t(nn,nnn)**2)nnn=nnnn 1040 continue ipt(nn)=nnn 1031 continue write(4,'(i4,'':'',6I12)') * (nn,(ipt(nnn),nnn=nn,nn+5),nn=1,mstot,6) ! dump c write(4,'(6I12)')(ipt(nn),nn=1,mstot) ! dump write(4,'(1x)') ! dump endif ! dump C if(mstot.gt.msort)then write(*,'(1x// * '' ERROR: the current MSORT='',i5,'' is too small''/ * '' and needs to be increased''//)') * msort stop endif DO 1020 I=1,MStot WK(I)=H1(I,I) IPT(I)=I 1020 CONTINUE c c...subsort eigenvalues within each J-block in each F-block in ascending c order of Tau=(K_a)-(K_c) assuming that diagonalisation c 1/ does not alter the K positions within a J block c 2/ does not introduce swaps between F blocks c c (the loop deals with multiple M.F blocks, although for one quadrupole c only one unconnected M.F block is set up at a time, ie mfmin=mfmax) c mm=0 i=I1+I2 do 3003 mfloop=mfmin,mfmax,2 ! Fortran2018 if(idump.ne.0)write(4,1030)mfloop,VCALC 1030 format(' ----------- M.F = ',i2,'/2 ------------',F11.4,' V:') do 2003 f=flow,fhigh,2 ! Fortran2018 jh=F+i jl=iabs(i-F) DO 1003 j=jl,jh,2 IF(IDUMP.gE.2) ! dump * write(4,1005)i,F,j/2,(h1(nn,nn),nn=mm+1,mm+(j+1)) ! dump CALL SORTJ(j,mm+1,mm+j+1) ! <---- SORTJ IF(IDUMP.ge.1) ! dump * write(4,1025)i,F,j/2,(h1(IPT(nn),IPT(nn)), ! dump * nn=mm+1,mm+j+1) ! dump mm=(j+1)+mm if(mm.ge.mstot)goto 1015 1003 continue 2003 continue ! Fortran2018 3003 continue ! Fortran2018 1005 format(/1X,'I=',I2,'/2 F=',I2,'/2 J=',I2,' RAW E-levels:', ! dump * 20(/1X,11F13.4)) ! dump 1025 format(1X,'I=',I2,'/2 F=',I2,'/2 J=',I2,' SORTED LEVELS:', ! dump * 20(/1X,11F13.4)) ! dump C C C...determine which eigenvalues will be necessary and load the relevant C eigenvalues into TFREQ(i,2) with the appropriate sign C 1015 NEXTR=0 DO 102 LL=istart,ifin NLEV=2 IF(INFOE(LL,1).LT.0)NLEV=1 NLINE=IABS(INFOE(LL,1)) I=LABO(NLINE,NLEV,1) J=LABO(NLINE,NLEV,2) K=LABO(NLINE,NLEV,3) F=labo(nline,nlev,4) c CALL NUMENL(I1,I2,F,I,J,K,MFCALC,MFMIN,MFMAX, ! <---- NUMENL * FLOW,FHIGH,neigv) M=IPT(neigv) c NEXTR=NEXTR+1 ! old label 120 NELINF(NEXTR,1)=INFOE(LL,1) NELINF(NEXTR,2)=M if(nline.gt.nt.or.nline.lt.1.or. * m.gt.mdim.or.m.lt.1.or.ll.gt.2*nt)then write(3,360)nline,m,ll,nextr,neigv,i1,i2,F,i,j,k 360 format(1x/' INDEX ERROR:'/20x,'nline=',i8/ * 20x,' m=',i8/20x,' ll=',i8/20x,'nextr=',i8/ * 20x,'neigv=',i8//' i1,i2,f,i,j,k: ',6i4//) stop endif TFREQ(NLINE,2)=TFREQ(NLINE,2)+H1(M,M)*ISIGN(1,INFOE(LL,1)) 102 CONTINUE C C...evaluate the derivatives for the levels extracted above (note that C in the least-squares procedure the sign of correction for a constant C which is in an imaginary term (chi.ac, chi.bc, mu.c) is reversed C DO 107 LL=1,NCON ! old label 106 DO 108 I=1,NCONST CONST(I)=0.0D0 108 CONTINUE lll=icfit(ll) CONST(LLL)=1.d0 if(lll.ge.32.and.lll.le.34)const(35)=con(35) if(lll.eq.20.or.lll.eq.21)then const(17)=1.D-15 endif if(lll.eq.22)then const(18)=1.D-15 endif if(lll.eq.35)then const(35)=-con(35)*con(35) do 444 i=32,34 const(i)=con(i) 444 continue endif C CALL SETHAM(MStot,I1,FLOW,FHIGH,CONST,vcalc,mfmin,mfmax) ! <---- SETHAM C DO 109 I=1,NEXTR NLINE=IABS(NELINF(I,1)) M=NELINF(I,2) EE=0.0 DO 110 M1=1,MStot DO 111 M2=M1,MStot RMULT=2.D0 IF(M1.EQ.M2)RMULT=1.D0 c if(imag.eq.1)then C C...(a+ib)(c+id) = ac-bd+i(bc+ad) C TR=T(M1,M)*T(M2,M)-T(M1+MStot,M)*T(M2+MStot,M) C TI=T(M1+MStot,M)*T(M2,M)+T(M1,M)*T(M2+MStot,M) C TR=H1(M1,M2)*TR+H1(M1+MStot,M2)*TI C...(a+ib)(c+id)* = (a+ib)(c-id)* = ac+bd+i(bc-ad) TR=T(M1,M)*T(M2,M)+T(M1+MStot,M)*T(M2+MStot,M) TI=T(M1+MStot,M)*T(M2,M)-T(M1,M)*T(M2+MStot,M) TR=H1(M1,M2)*TR-H1(M1+MStot,M2)*TI EE=EE+RMULT*TR else EE=EE+RMULT*H1(M1,M2)*T(M1,M)*T(M2,M) endif c 111 CONTINUE 110 CONTINUE DER(NLINE,LL)=DER(NLINE,LL)+ISIGN(1,NELINF(I,1))*EE 109 CONTINUE 107 CONTINUE C istart=ifin+1 ! old label 115 if(istart.le.nelev)goto 502 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c RETURN END c c_____________________________________________________________________________ c FUNCTION MHIGH(I1,I2,F) C C MHIGH returns the dimension of the eigenvalue matrix H for given I1,I2 C and F where: C C I1+I2 I+F C \--- \--- C MHIGH = > > (2J+1) C /--- /--- C I=|I1-I2| J=|I-F| C C eg. - for I1=1, I2=3/2 and Fmax=13/2 there will be: C C I = 1/2: J = 6,7 C I = 3/2: J = 5,6,7,8 C I = 5/2: J = 4,5,6,7,8,9 ---> MHIGH = 168 C C - for I1=3/2, I2=0, Fmax=41/2 ---> MHIGH = 168 C - for I1=3/2, I2=0, Fmax=165/2 (J=80) ---> MHIGH = 664 C........................................................................ C Convenient program size estimator: C C PROGRAM SIZE (bytes) = nmax*nmax*12 + (200000 to 400000) C C where: nmax=(2I1+1)(2I2+1)(2F+1), and the 0.2-0.4Mb element C depends on system and FORTRAN used C........................................................................ C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) c MS=0 IH=I1+I2 IL=IABS(I1-I2) DO 1000 I=IL,IH,2 ! Fortran2018 JH=F+I JL=IABS(I-F) DO 1 J=JL,JH,2 MS=(J+1)+MS 1 CONTINUE 1000 CONTINUE ! Fortran2018 MHIGH=MS c RETURN END c c_____________________________________________________________________________ c SUBROUTINE QUANTN(M,I,F,J,K,MF,MFMIN,MFMAX,FLOW,FHIGH) C C Called from SETHAM() C C This routine returns F,J,K,MF for the M'th row or column of the eigenvalue C matrix set up for total nuclear spin I and containing C 1. M.F blocks from MFMIN to MFMAX C 2. in each M.F block F blocks from F=MF to FHIGH C C NOTE: I, F, J, MF, MFMIN, MFMAX, FHIGH are all in units of 1/2 C K is in units of 1 C C........................................................................ c c Blocks in the J,K,I,F,M.F basis: c c I=1 c c F M.F J 2J+1 Sum SumSum c c 0 0 1 3 3 3 c 1 0,1 0,1,2 1, 3, 5 9 12 c 2 0,1,2 1,2,3 3, 5, 7 15 27 c 3 0,1,2,3 2,3,4 5, 7, 9 21 49 c 4 0,1,2,3,4 3,4,5 7, 9,11 27 76 c 5 0,1,2,3,4,5 4,5,6 9,11,13 33 109 c 6 0,1,2,3,4,5,6 5,6,7 11,13,15 39 148 c 7 0,1,2,3,4,5,6,7 6,7,8 13,15,17 45 193 c 8 0,1,2,3,4,5,6,7,8 7,8,9 15,17,19 51 244 C c I=3/2 c c 2F 2M.f J 2J+1 Sum SumSum c c 1 1 1,2 3, 5 8 8 c 3 1,3 0,1,2,3 1, 3, 5, 7 16 24 c 5 1,3,5 1,2,3,4 3, 5, 7, 9 24 48 c 7 1,3,5,7 2,3,4,5 5, 7, 9,11 32 80 c 9 1,3,5,7,9 3,4,5,6 7, 9,11,13 40 120 c 11 1,3,5,7,9,11 4,5,6,7 9,11,13,15 48 168 c 13 1,3,5,7,9,11,13 5,6,7,8 11,13,15,17 56 224 c 15 1,3,5,7,9,11,13,15 6,7,8,9 13,15,17,19 64 288 c 17 1,3,5,7,9,11,13,15,17 7,..,10 15,17,19,21 72 360 c 19 1,................,19 8,..,11 17,19,21,23 80 440 c 21 1,................,21 9,..,12 19,21,23,25 88 528 c 23 1,................,23 10,.,13 21,23,25,27 96 624 c 25 1,................,25 11,.,14 23,25,27,29 104 728 C c I=5/2 c c 2F 2M.f J 2J+1 Sum SumSum c c 1 1 2,3 5, 7 12 12 c 3 1,3 1,2,3,4 3, 5, 7, 9 24 36 c 5 1,3,5 0,1,2,3,4,5 1, 3, 5, 7, 9,11 36 72 c 7 1,3,5,7 1,2,3,4,5,6 3, 5, 7, 9,11,13 48 120 c 9 1,3,5,7,9 2,3,4,5,6,7 5, 7, 9,11,13,15 60 180 c 11 1,....,11 3,4,5,6,7,8 7, 9,11,13,15,17 72 252 c 13 1,....,13 4,5,6,7,8,9 9,11,13,15,17,19 84 336 c 15 1,....,15 5,......,10 11,13,15,17,19,21 96 432 c 17 1,....,17 6,......,11 13,15,17,19,21,23 108 540 c 19 1,....,19 7,......,12 15,17,19,21,23,25 120 660 c 21 1,....,21 8,......,13 17,19,21,23,25,27 132 792 c 23 1,....,23 9,......,14 19,21,23,25,27,29 144 936 c 25 1,....,25 10,......,15 21,23,25,27,29,31 156 1092 c 27 1,....,27 11,......,16 23,25,27,29,31,33 168 1260 C 29 1,....,29 12,......,17 25,27,29,31,33,35 180 1440 C 31 1,....,31 13,......,19 27,29,31,33,35,37 192 1632 C 33 1,....,33 15,......,21 29,31,33,35,37,39 204 1836 C C........................................................................ C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) c M1=M DO 2000 MFLOOP=MFMIN,MFMAX,2 ! Fortran2018 DO 1000 F=FLOW,FHIGH,2 JH=F+I JL=IABS(F-I) DO 1 J=JL,JH,2 M2=M1-J-1 IF(M2.LE.0)GOTO 2 M1=M2 1 CONTINUE 1000 CONTINUE ! Fortran2018 2000 CONTINUE ! Fortran2018 c 2 K=M1-1-J/2 MF=MFLOOP C RETURN END c c_____________________________________________________________________________ c SUBROUTINE NUMENL(I1,I2,F,I,J,K,MF,MFMIN,MFMAX,FLOW,FHIGH,M) C C Called from READ() C C This routine returns M, the number of eigenvalue among the eigenvalues of C H, for given I,J,K=(Tau+J+1) and I1,I2,F,MF containing M.F blocks C from MFMIN to MFMAX, each for F from F=M.F to F=FHIGH C C NOTE: I1,I2,F,FHIGH,MFMIN,MXMAX and I are all in units of 1/2, whereas C J and K are in units of 1 C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) C J2=J*2 I=I1+I2 M=0 do 1000 MFLOOP=MFMIN,MFMAX,2 ! Fortran2018 DO 1 FF=FLOW,FHIGH,2 JH=FF+I JL=IABS(FF-I) DO 2 JJ1=JL,JH,2 IF(MFLOOP.LT.MF)GOTO 3 IF(FF.LT.F)GOTO 3 IF(JJ1.LT.J2)GOTO 3 M=M+K GOTO 5 3 M=M+JJ1+1 2 CONTINUE 1 CONTINUE 1000 CONTINUE ! Fortran2018 C 5 RETURN END c c_____________________________________________________________________________ c c SUBROUTINE SETHAM(MStot,I1,FLOW,FHIGH,CON,vcalc,mfmin,mfmax) C C Called from READ() to set up the H matrix of size MStot x MStot: C C - for nucleus with spin I1 C - containing all M.F blocks from MFMIN to MFMAX C - for each M.F block containing all F blocks from the lowest value C F=M.F up to F=FHIGH C - from constants in CON c - for Stark voltage VCALC C C Chi.ac, Chi.bc and mu.c give rise to imaginary components in the c Hamiltonian. If any of these are the only off-diagonal components c of this type (i.e. totally imaginary off-diagonals, then it is c possible to use a standard calculation - the known effect of Chi.ac c being fittable as a chi.ab constant). If off-diagonals contain c a mixture of real and imaginary elements then it is necessary to switch c to the full imaginary treatment. C C In that case the Hermitian H matrix is diagonalised via the augmented C real approach ie. the complex elements (a+ib) are expressed using C the real matrix C H.aug = A -B C B A C Diagonalisation of H.aug yields two sets of identical eigenvalues C and eigenvectors such that first MStot elements in a given column are C the real components, with the following MStot elements being the imaginary C components. C use globdat C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=36,NLIN=2000) COMMON /HAMBLK/D,G,GI,GJ,GK,V,SR1,SR2,dip,dipi,dspac,CONST common /stark/volts(nlin) DIMENSION CON(NCONST),CONST(NCONST) DIMENSION D(5),G(5),GI(5),GJ(5),GK(5),V(5),SR1(5),SR2(5), * dip(5),dipi(5) ROOT6=DSQRT(6.D0) ROOT23=DSQRT(2.D0/3.D0) root12=dsqrt(1.d0/2.d0) C C...Use the constants to set up the following tensors: C C D(5) - spherical spin-spin tensor C G(5) - spherical quadrupole tensor for nucleus with spin I1, real part C GI(5) - imaginary part of the spherical quadrupole tensor (for chi.bc C and chi.ac) C GJ(5) - J(J+1) term on G C GK(6) - K**2 term on G C V(5) - spherical quadrupole tensor for nucleus with spin I2 C SR1(5) - spherical spin-rotation tensor for nucleus with spin I1 C SR2(5) - spherical spin rotation tensor for nucleus with spin I2 C dip(5) - spherical dipole moment tensor (this only has three non-zero c elements -1,0,1) - real part only C dipi(5) - imaginary part of the spherical dipole moment tensor, the only C non-zero elements are -1 and 1 from mu.c C C Order of components in all tensors: C(-2),C(-1),C(O),C(1),C(2) C DO 30 I=1,5 D(I)=0.0D0 G(I)=0.0D0 GI(I)=0.d0 GJ(I)=0.0D0 GK(I)=0.0D0 V(I)=0.0D0 SR1(I)=0.0D0 SR2(I)=0.0D0 dip(i)=0.d0 dipi(i)=0.d0 30 CONTINUE ! chi_spherical: C ! G(1)= CON(18)/ROOT6 ! -2 G(2)= ROOT23*CON(19) ! chi_ab (xz) -1 G(3)=CON(17) ! 0 G(4)= -ROOT23*CON(19) ! chi_ab 1 G(5)= CON(18)/ROOT6 ! 2 c ! GI(1)=-2.D0*CON(27)/ROOT6 ! chi_bc (xy) -2 GI(2)=-ROOT23*CON(28) ! chi_ac (yz) -1 GI(3)= 0.d0 ! 0 GI(4)=-ROOT23*CON(28) ! chi_ac 1 GI(5)= 2.D0*CON(27)/ROOT6 ! chi_bc 2 C GJ(3)=CON(20) GJ(1)=CON(22)/ROOT6 GJ(5)=CON(22)/ROOT6 C GK(3)=CON(21) C c V(3)=CON(20) c V(1)=CON(21)/ROOT6 c V(5)=CON(21)/ROOT6 c V(2)=ROOT23*CON(22) c V(4)=-ROOT23*CON(22) C SR1(3)=CON(23)-CON(25)/3.D0 SR1(1)=CON(24)/ROOT6 SR1(5)=CON(24)/ROOT6 C SR2(3)=CON(26)-CON(28)/3.D0 SR2(1)=CON(27)/ROOT6 SR2(5)=CON(27)/ROOT6 C D(3)=CON(29) D(1)=CON(30)/ROOT6 D(5)=CON(30)/ROOT6 D(2)=ROOT23*CON(31) D(4)=-ROOT23*CON(31) C dip(3)=con(32) c dip(2)= root12*(con(33)-con(34)) c dip(4)=-root12*(con(33)+con(34)) c dip(2)= root12* con(33) dip(4)= -root12* con(33) dipi(2)= root12*( -con(34)) dipi(4)=-root12*( +con(34)) c dspac=con(35) c DO 102 L=1,NCONST CONST(L)=CON(L) 102 CONTINUE C C...Set up H, assuming it is Hermitian C DO 1101 L=1,MStot ! Fortran2018 DO 101 M=L,MStot CALL QUANTN(L,I1,FR,JR,KR,MFR,mfmin,mfmax,FLOW,fhigh) ! <---- QUANTN CALL QUANTN(M,I1,FC,JC,KC,MFC,mfmin,mfmax,FLOW,fhigh) ! <---- QUANTN C write(4,'('' row='',i3,'' F,J,K,MF'',4i3,3x, C * '' col='',i3,'' F,j,k,MF'',4i3)') C * L,fr,jr,kr,mfr, M,fc,jc,kc,mfc CALL HAMIL(HREAL,HIMAG,I1,I2,FR,FC,JR,JC,KR,KC,MFR,MFC,vcalc) ! <---- HAMIL H1(L ,M) =HREAL H1(L+MStot,M+MStot)=HREAL H1(l ,m+mstot)= himag h1(l+mstot,m) =-himag IF(L.NE.M)THEN H1(M ,L) =HREAL H1(M+MStot,L+MStot)=HREAL H1(M ,L+MStot)=-HIMAG H1(M+MStot,L) = HIMAG ENDIF C IF(L.NE.M)H1(M,L)=H1(L,M) 101 CONTINUE 1101 CONTINUE ! Fortran2018 C RETURN END C C_____________________________________________________________________________ C SUBROUTINE HAMIL(HREAL,HIMAG, * I1,I2,FR,FC,JR,JC,KR,KC,MFR,MFC,vcalc) C C Called: M/P -> READ() -> SETHAM() -> HAMIL() C C Evaluation of a Hamiltonian element for given I1,I2,F and: C C row F,I,J,K,MF quantum numbers: FR,JR,KR,MFR C column F,I,J,K,MF quantum numbers: FC,JC,KC,MFC C C I1,I2,FR,IF,JR,JC,MFR,MFC are all in units of 1/2 C KR,KC are in units of 1 C C quantum numbers fed into the THREEJ() etc. routines C are all in units of 1: FF, MF, I, J, K C FFPR,MFPR,IPR,JPR, KPR C and they are all real quantities C C use globdat C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=36,NLIN=2000) COMMON /HAMBLK/D,G,GI,GJ,GK,H,SR1,SR2,dip,dipi,dspac,CON common /stark/volts(nlin) DIMENSION D(5),G(5),GI(5),GJ(5),GK(5),H(5),SR1(5),SR2(5), * dip(5),dipi(5),CON(NCONST) REAL(8) J,JPR,K,KPR,IA,IB,Q,I,IPR,FF,FFPR,JJP1,KSQ,mf,mfpr c S=0.0D0 SS=0.0D0 ZERO=0.d0 ONE=1.d0 TWO=2.d0 IR=I1 IC=I1 C IF(MFR.NE.MFC.OR.FR.NE.FC.OR.JR.NE.JC)GOTO 1 F=FR JR=JR/2 JJP1=JR*(JR+1) KSQ=KR*KR IF(KR.EQ.KC)GOTO 2 IF(KR.EQ.KC+2.OR.KR.EQ.KC-2)GOTO 3 GOTO 7 C C C 2 HBPC=0.5D0*(CON(2)+CON(3)) if(LINSYM.eq.1)hbpc=2.d0*hbpc S=S+HBPC*JJP1+(CON(1)-HBPC)*KSQ-CON(4)*(JJP1*JJP1) 1 -CON(5)*(JJP1*KSQ)-CON(6)*(KSQ*KSQ) 1 +CON(9)*JJP1**3+CON(10)*JJP1*JJP1*KSQ 1 +CON(11)*(JJP1*KSQ*KSQ)+CON(12)*KSQ*KSQ*KSQ c c Chi.d term according to the symmetric top definition c if(con(26).ne.0.d0.and.JR.ne.0.d0)then I=DBLE(IR)*0.5D0 FF=DBLE(FR)*0.5D0 C=FF*(FF+1.d0)-I*(I+1.d0)-JJP1 YIJF=(0.75D0*C*(C+1.d0)-I*(I+1.D0)*JJP1)/ 1 (2.d0*(2.d0*JR-1.d0)*(2.d0*JR+3.d0)*I*(2.d0*I-1.d0)) S=S+(KSQ*(4.d0*KSQ-1.d0)/JJP1)*YIJF*CON(26) endif c GOTO 7 C C C 3 HBMC=CON(2)-CON(3) if(LINSYM.EQ.1)HBMC=0.0D0 if(kr.eq.kc-2)then S=S+(0.25D0*HBMC-CON(7)*JJP1-0.5D0*CON(8) 1 *((KR+2)**2+KSQ)+CON(13)*JJP1*JJP1 1 +0.5D0*CON(14)*JJP1*((KR+2)**2+KSQ) 1 +0.5D0*CON(15)*((KR+2)**4+KSQ*KSQ)) 2 *DSQRT(DBLE((JJP1-KR*(KR+1))*(JJP1-(KR+1)*(KR+2)))) else S=S+(0.25D0*HBMC-CON(7)*JJP1-0.5D0*CON(8) 1 *((KR-2)**2+KSQ)+CON(13)*JJP1*JJP1 1 +0.5D0*CON(14)*JJP1*((KR-2)**2+KSQ) 1 +0.5D0*CON(15)*((KR-2)**4+KSQ*KSQ)) 2 *DSQRT(DBLE((JJP1-KR*(KR-1))*(JJP1-(KR-1)*(KR-2)))) endif C C 7 JR=2*JR C 1 J=DBLE(JR)*0.5D0 JPR=DBLE(JC)*0.5D0 K=KR KPR=KC IA=DBLE(I1)*0.5D0 IB=DBLE(I2)*0.5D0 I=DBLE(IR)*0.5D0 IPR=DBLE(IC)*0.5D0 FF=DBLE(FR)*0.5D0 FFPR=dble(FC)*0.5d0 c Q=kc-kr M=-nint(Q) mf=dble(mfr)*0.5d0 mfpr=dble(mfc)*0.5d0 C C C C Matrix element is -Field*Eq.(14) from ref.(H.E). Since only the mu.Zg C terms are required the value of r in the first threej symbol is set C to zero. C C...The conversion factor for (E*mu) of 0.50341125 is to use C field in V/cm, dipoles in Debyes and frequency in MHz: C 0.50341125 = 10/(2.99792458*6.6260755); 1986 constants, C and only the numerical value of the Planck constant is subject to change C C Note that the Debye was defined in the CGS era as 10^(-10) of the C dipole moment from two charges equal to +- 1 cgs charge unit (franklin) C separated by 1 Angstrom, thus 1D = 10^(-18) franklin cm; C 1 Franklin = (1/c(CGS)) Coulomb = 3.33564095*10^(-10) C so that C 1D = 3.33564095*10^(-30) C m C if(dspac.ne.0.0d0.and.iabs(nint(q)).le.1)then C n=(fr+fc+mfr+ir)/2+kr+1 n=nint(ff+ffpr+mf+i+k+one) TERM=-(-1)**(n) * dsqrt(dble( (jr+1)*(jc+1)*(fr+1)*(fc+1) )) 1 *sixj ( j, ff, i, 1 ffpr, jpr, one ) 1 *threej( ff, one, ffpr, 1 -mf, zero, mfpr) 1 *threej( j, one, jpr, 1 -k, -q , kpr ) c S=S+TERM* 0.50341125d0*vcalc*(1.d0+con(36)*vcalc)* 1 dip(3+NINT(q))/dspac SS=SS+TERM*0.50341125d0*vcalc*(1.d0+con(36)*vcalc)* 1 dipi(3+NINT(q))/dspac endif c if(mfr.ne.mfc.or.fr.ne.fc)goto 9 C C M=M+3 IF(M.GT.5.OR.M.LT.1)then S=0.0D0 GOTO 9 ENDIF C C C C Note that: C C SIXJ| I, IPR, TWO | * SIXJ| I, IPR,TWO| = [(-1)^(2+I+I)/(2I+1)] C | JPR, J, FF| | IA, IA, IB| C C * SIXJ| FF, I, J| C |TWO, JPR, I| C C so that the double nucleus formula below reduces to eq.(9) of ref.(H.E) C IF(I1.LE.1)GOTO 5 IF(G(M).EQ.0.0D0.AND.GI(M).EQ.0.D0)GOTO 5 C N=IC+(F+I1+I2)/2+KR TERM=(-1)**N*0.25D0*DSQRT(DBLE((IR+1)*(IC+1)*(JR+1)*(JC+1))) 1 * SIXJ( I, IPR, TWO, 1 JPR, J, FF) 2 *THREEJ( J,TWO, JPR, 2 -K, -Q, KPR ) 3 *( SIXJ( I, IPR,TWO, 3 IA, IA, IB) 4 /THREEJ( IA, TWO, IA, 4 -IA,ZERO, IA ) 5 ) c S= S+TERM*( G(M) + GJ(M)*J*(J+1) + GK(M)*KR*KR ) SS=SS+TERM* GI(M) C C C 5 IF(I2.LE.1)GOTO 10 IF(H(M).EQ.0.0D0)GOTO 10 C N=KR+(F+IC+IR+I1+I2)/2 S=S+(-1)**N*0.25D0*DSQRT(DBLE((IR+1)*(IC+1)*(JR+1)*(JC+1))) 1 * SIXJ( I, IPR, TWO, 1 JPR, J, FF) 2 *THREEJ( J,TWO, JPR, 2 -K, -Q, KPR ) 3 *( SIXJ( I, IPR,TWO, 3 IB, IB, IA) 4 /THREEJ( IB, TWO, IB, 4 -IB,ZERO, IB ) ) 5 *H(M) C C C 10 IF(D(M).EQ.0.0D0)GOTO 6 C N=KR+(IC+F)/2 S=S+(-1)**N*DSQRT((2.0D0*I+1.0D0)*(2.0D0*IPR+1.0D0)*(2.0D0*J 1+1.0D0)*30.D0*(2.0D0*JPR+1.0D0)*(2.0D0*IA+1.0D0)*(IA+1.0D0)*IA 2*(2.0D0*IB+1.0D0)*(IB+1.0)*IB) 3 *THREEJ( J,TWO,JPR, -K,-Q,KPR) 4 * SIXJ( I,IPR,TWO, JPR, J, FF) 5 * X9J(IA, IA,ONE, IB,IB,ONE, I,IPR,TWO)*D(M)/2.0D0 C C C 6 IF(SR1(M).EQ.0.0D0)GOTO 11 C N=IC+KR+(I1+I2+JR+JC+F)/2 S=S+(-1)**N*DSQRT(30.0D0*(2.0D0*J+1.0D0)*(2.0D0*JPR+1.0D0) 1*(2.0D0*J+1.0D0)*(J+1.0D0)*J*(2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 2*(2.0D0*IA+1.0D0)*(IA+1.0D0)*IA)*SIXJ(IA,I,IB,IPR,IA,ONE) 3 * SIXJ( I, J, FF, JPR,IPR,ONE) 4 * SIXJ(ONE,TWO,ONE, JPR, J, J) 5 *THREEJ( J,TWO,JPR, -K, -Q,KPR)*SR1(M)/2.0D0 C 11 IF(JR.NE.JC.OR.KR.NE.KC)GOTO 8 ! C.D.Esp IF(CON(25).EQ.0.0D0)GOTO 8 C N=(I1+I2+F+JR)/2+IC+1 S=S-(-1)**N*DSQRT((2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 1*(2.0D0*IA+1.0)*(IA+1.0D0)*IA*(2.0D0*J+1.0D0)*(J+1.0D0)*J) 2 *SIXJ( I,J,FF, J,IPR,ONE) 3 *SIXJ(IA,I,IB, IPR, IA,ONE) 4 *CON(25)/3.D0 ! C.D.Esp C C C 8 IF(SR2(M).EQ.0.0D0)GOTO 12 C N=(IC+JR+F+I1+I2+IR+JC)/2+KR S=S+(-1)**N*DSQRT(30.D0*(2.0D0*J+1.0D0)*(2.0D0*JPR+1.0D0) 1 *(2.0D0*J+1.0D0)*(J+1.0D0)*J*(2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 2 *(2.0D0*IB+1.0D0)*(IB+1.0D0)*IB)*SIXJ(I,ONE,IPR,IB,IA,IB) 3 * SIXJ( I, J, FF, JPR,IPR,ONE) 4 * SIXJ(ONE,TWO,ONE, JPR, J, J) 5 *THREEJ( J,TWO,JPR, -K, -Q,KPR)*SR2(M)/2.0D0 C 12 IF(JR.NE.JC.OR.KR.NE.KC)GOTO 9 ! C.D.Esp IF(CON(28).EQ.0.0D0)GOTO 9 C N=(I1+I2+IR+IC+JR+F)/2+1 S=S-(-1)**N*DSQRT((2.0D0*IPR+1.0D0)*(2.0D0*I+1.0D0) 1 *(2.0D0*IB+1.0D0)*(IB+1.0D0)*IB*(2.0D0*J+1.0D0)*(J+1.0D0)*J) 2 *SIXJ(I, J, FF, J,IPR,ONE) 4 *SIXJ(I,ONE,IPR, IB, IA, IB) 5 *CON(28)/3.D0 ! C.D.Esp C C 9 IF(ABS(S). LE.1.D-300)S=0.D0 IF(ABS(SS).LE.1.D-300)SS=0.D0 HREAL=S HIMAG=SS c RETURN END c c_____________________________________________________________________________ c C C ANGULAR MOMENTUM MATRIX ELEMENT CALCULATIONS C c c_____________________________________________________________________________ c c FUNCTION THREEJ(X1,X2,X3,Y1,Y2,Y3) IMPLICIT REAL(8) (A-H,O-Z) c IF(ABS(Y1+Y2+Y3).GT.0.1D0)GOTO 50 TR=TRIANG(X1,X2,X3) IF(NINT(TR).NE.0)GOTO 60 50 THREEJ=0.0D0 itemp=int( abs(10.*(Y1+Y2+Y3)) ) if(itemp.le.6.and.itemp.gt.1)write(*,'(1x,a,3F10.3)') * ' THREEJ warning:',y1,y2,y3 RETURN 60 IR=int( X2+Y3-X1 ) PP=(-1)**IR QQ=TRI(X1,X2,X3) RR= DSQRT(FAC(X1-Y1))*DSQRT(FAC(X1+Y1))*DSQRT(FAC(X2-Y2)) RR=RR*DSQRT(FAC(X2+Y2))*DSQRT(FAC(X3-Y3))*DSQRT(FAC(X3+Y3)) SS=SUML(X1,X2,X3,Y1,Y2,Y3) THREEJ=QQ*RR*SS*PP c RETURN END c c_____________________________________________________________________________ c FUNCTION SIXJ(X1,X2,X3,Y1,Y2,Y3) IMPLICIT REAL(8) (A-H,O-Z) c SIXK=TRIANG(X1,X2,X3)*TRIANG(Y1,X2,Y3)*TRIANG(X1,Y2,Y3)* * TRIANG(Y1,Y2,X3) IF(NINT(SIXK).EQ.0)GOTO 2 IR=int( X1+X2+Y1+Y2 ) PP=(-1)**IR SIXK=PP*TRI(X1,X2,X3)*TRI(Y1,Y2,X3)*TRI(Y1,X2,Y3)*TRI(X1,Y2,Y3) 1 *SUMK(X1,X2,X3,Y1,Y2,Y3) 2 CONTINUE SIXJ=SIXK c RETURN END c c_____________________________________________________________________________ c FUNCTION X9J(X1,X2,X3,X4,X5,X6,X7,X8,X9) IMPLICIT REAL(8) (A-H,O-Z) c S=0.0D0 Y1=TRIANG(X1,X2,X3) Y2=TRIANG(X4,X5,X6) Y3=TRIANG(X7,X8,X9) Y4=TRIANG(X1,X4,X7) Y5=TRIANG(X2,X5,X8) Y6=TRIANG(X3,X6,X9) TT=Y1*Y2*Y3*Y4*Y5*Y6 IF(TT.LT.0.1D0)GOTO 1 Z1=X1+X9 Z2=X6+X2 Z3=X4+X8 N1=MINT(Z1) N2=MINT(Z2) N3=MINT(Z3) IF(N1.NE.N2)GOTO 1 IF(N2.NE.N3)GOTO 1 IF(N1.EQ.0)GOTO 2 DO 100 I=0,10 X0=I S=S+(2.0D0*X0+1.0D0)*SIXJ(X1,X2,X3,X6,X9,X0) 1 *SIXJ(X4,X5,X6,X2,X0,X8)*SIXJ(X7,X8,X9,X0,X1,X4) 100 CONTINUE GOTO 1 2 DO 101 I=1,19,2 X0=I/2.0D0 S=S-(2.0D0*X0+1.0D0)*SIXJ(X1,X2,X3,X6,X9,X0) 1 *SIXJ(X4,X5,X6,X2,X0,X8)*SIXJ(X7,X8,X9,X0,X1,X4) 101 CONTINUE 1 X9J=S c RETURN END c c_____________________________________________________________________________ c FUNCTION TRIANG(O,P,Q) C C Test for the triangular condition: this is not satisfied if sum of any C two of the three numbers is smaller than the third number C IMPLICIT REAL(8) (A-H,O-Z) c Z0=0.d0 IF((O+P-Q).LT.Z0.OR.(O-P+Q).LT.Z0.OR.(-O+P+Q).LT.Z0)GOTO 1 GOTO 2 1 STRIA=Z0 TRIANG=STRIA RETURN 2 STRIA=1.0D0 TRIANG=STRIA c RETURN END c c_____________________________________________________________________________ c FUNCTION TRI(U,V,W) IMPLICIT REAL(8) (A-H,O-Z) C Z1=1.0D0 TRK=DSQRT((FAC(U+V-W)*FAC(U-V+W)*FAC(-U+V+W))/(FAC(U+V+W+Z1))) TRI=TRK C RETURN END c c_____________________________________________________________________________ c FUNCTION FAC(Z) IMPLICIT REAL(8) (A-H,O-Z) c X=Z Z1=1.0D0 FAC=Z1 1 IF(X.LE.Z1)GOTO 3 FAC=FAC*X ! old label 2 X=X-Z1 GOTO 1 3 CONTINUE c RETURN END c c_____________________________________________________________________________ c FUNCTION SUML(U,V,W,X,Y,Z) IMPLICIT REAL(8) (A-H,O-Z) c DIMENSION FF(5) FF(1)=U-X FF(2)=V+Y FF(3)=U+V-W FF(4)=W-V+X FF(5)=W-U-Y dummy=z DO 75 IK=2,3 IF(FF(IK).GE.FF(1))GOTO 75 FF(1)=FF(IK) 75 CONTINUE KMAX=int( FF(1) ) IF(FF(4).LE.FF(5))GOTO 80 FF(4)=FF(5) 80 IF(FF(4).LE.0.0D0)GOTO 85 KMIN=0 GOTO 95 85 KMIN=int( -FF(4) ) 95 SUMQ=0.0D0 DO 97 KF=KMIN,KMAX XY=(-1)**KF RK=KF SUMQ=SUMQ+XY/(FAC(U-X-RK)*FAC(W-V+X+RK)*FAC(V+Y-RK)* 1 FAC(W-U-Y+RK)*FAC(RK)*FAC(U+V-W-RK)) 97 CONTINUE SUML=SUMQ c RETURN END c c_____________________________________________________________________________ c FUNCTION SUMK(S1,S2,S3,T1,T2,T3) IMPLICIT REAL(8) (A-H,O-Z) DIMENSION FA(7) c Z1=1.0D0 FA(1)=S1+S2+T1+T2+Z1 FA(2)=S1+S2-S3 FA(3)=T1+T2-S3 FA(4)=S1+T2-T3 FA(5)=T1+S2-T3 FA(6)=-S1-T1+S3+T3 FA(7)=-S2-T2+S3+T3 DO 1 IK=2,5 IF(FA(IK).LT.FA(1))GOTO 2 GOTO 1 2 FA(1)=FA(IK) 1 CONTINUE KFMAX=int( FA(1) ) IF(FA(6).LE.FA(7))GOTO 3 FA(6)=FA(7) ! old label 4 Z0=0.0D0 3 IF(FA(6).LE.Z0)GOTO 5 GOTO 6 5 KFMIN=int( -FA(6) ) GOTO 8 6 KFMIN=0 8 SUMP=Z0 DO 7 KF=KFMIN,KFMAX XY=(-1)**KF RK=KF Z1=1.0D0 SUMP=SUMP+XY*FAC(S1+S2+T1+T2+Z1-RK)/(FAC(S1+S2-S3-RK)*FAC(T1+T2 1 -S3-RK)*FAC(S1+T2-T3-RK)*FAC(T1+S2-T3-RK)*FAC(-S1-T1+S3+T3+RK) 1 *FAC(-S2-T2+S3+T3+RK)*FAC(RK)) 7 CONTINUE SUMK=SUMP c RETURN END c c_____________________________________________________________________________ c FUNCTION MINT(A) IMPLICIT REAL(8) (A-H,O-Z) C C This is apparently a test for whether A is close to an integer (MINT=1) C or half integer (MINT=0) C I=int( 2.00001*A ) N=(1+(-1)**I)/2 MINT=N C RETURN END C c_____________________________________________________________________________ C C c_____________________________________________________________________________ c SUBROUTINE LSQRW(IN,LX,SD,X) C C Weighted least squares routine adapted from the LSQR routine from C A.S.Georgiou's program CDFIT C C The NINFIT variable was added to eliminate the problem C of counting unfitted lines in calculating the standard deviation) C C IN - number of equations C LX - number of parameters C SD - standard deviations on the fitted parameters C X - calculated increments to parameters C R - correlation matrix C WT - weight matrix, this is a vector containing the diagonal elements C of the M^-1 matrix of [ASZ] but without the common factor of C sigma_n^2, which cancels in the key Eq.28. In notes on the margin C the notional square form of this matrix is denoted by W C LINFIT - flag defining whether the line is (=1) or not (=0) to be used C in the fit C IMPLICIT REAL(8) (A-E,G-H,O-Z),INTEGER(F) PARAMETER (NCONST=36,NLIN=2000,NLIN2=2*NLIN) COMMON /BLKEN/TFREQ,DER,LABO,INFOE,volinf,ICFIT,NELEV COMMON /BLKINV/A COMMON /BLKFIT/WT,LINFIT,IFREE(nlin),R DIMENSION TFREQ(NLIN,3),DER(NLIN,NCONST),LABO(NLIN,2,5), * ICFIT(NCONST),INFOE(NLIN2,2),LINFIT(NLIN),WT(NLIN) DIMENSION A(NCONST,NCONST),R(NCONST,NCONST),C(NCONST), * SD(NCONST),X(NCONST),BC(NLIN),Z(NCONST),volinf(nlin2) c c...Deviation of fit: SIGMA = weighted [ASZ] Eq.30 sigma_w^2 c SIGMAF = unweighted c SIGMA=0.0D0 SIGMAF=0.d0 ninfit=0 DO 110 I=1,IN BC(I)=TFREQ(I,3) IF(LINFIT(I).EQ.1)then ninfit=ninfit+1 SIGMA=SIGMA+BC(I)*BC(I)*WT(I) SIGMAF=SIGMAF+BC(I)*BC(I) ENDIF 110 CONTINUE if(ninfit-LX.gt.0)then SIGMA=DSQRT(SIGMA/(ninfit-LX)) SIGMAF=DSQRT(SIGMAF/(ninfit-LX)) else SIGMA=DSQRT(SIGMA) SIGMAF=DSQRT(SIGMAF) endif C IFAIL=0 DO 1101 I=1,NCONST ! Fortran2018 DO 101 J=1,NCONST R(I,J)=0.0D0 101 CONTINUE 1101 CONTINUE ! Fortran2018 C C...Form: A = DER' * DER X'W X C DO 50 I1=1,LX DO 3 I2=I1,LX SUM=0.0D0 DO 5 IS=1,IN IF(LINFIT(IS).NE.1)GOTO 5 SUM=SUM+DER(IS,I1)*DER(IS,I2)*WT(IS) 5 CONTINUE A(I1,I2)=SUM 3 CONTINUE Z(I1)=DSQRT(1.D0/A(I1,I1)) 50 CONTINUE C C...rescale A to improve conditioning C DO 1051 I1=1,LX ! Fortran2018 DO 51 I2=I1,LX A(I1,I2)=A(I1,I2)*Z(I1)*Z(I2) A(I2,I1)=A(I1,I2) 51 CONTINUE 1051 CONTINUE ! Fortran2018 C C...Form: C = DER' * BC X'W Y C DO 6 I1=1,LX SUM=0.0D0 DO 8 IS=1,IN IF(LINFIT(IS).NE.1)GOTO 8 SUM=SUM+DER(IS,I1)*BC(IS)*WT(IS) 8 CONTINUE C(I1)=SUM 6 CONTINUE C C...Invert A ie. (DER'*DER) and 'unscale' (X'W X)^-1 C CALL INV(LX) ! <---- INV C DO 1052 I1=1,LX ! Fortran2018 DO 52 I2=1,LX A(I1,I2)=A(I1,I2)*Z(I1)*Z(I2) 52 CONTINUE 1052 CONTINUE ! Fortran2018 C C...Increments to constants X = A * C [ASZ],Eq.28: (X'W X)^-1 X W'Y C standard errors SD = SIGMA * A**1/2 C DO 9 I1=1,LX SUM=0.0D0 DO 10 IS=1,LX SUM=SUM+A(I1,IS)*C(IS) 10 CONTINUE X(I1)=SUM SD(I1)=SIGMA*DSQRT(A(I1,I1)) 9 CONTINUE C C...Correlation matrix R C SDEVSQ=SIGMA*SIGMA DO 1015 I1=1,LX ! Fortran2018 DO 15 I2=1,I1 R(I1,I2)=SDEVSQ*A(I1,I2)/(SD(I1)*SD(I2)) ! T_ij/(T_ii T_jj)^(1/2) 15 CONTINUE ! where T refers to Theta of [ASZ] 1015 CONTINUE ! Fortran2018 c RETURN END c c_____________________________________________________________________________ c SUBROUTINE HDIAGM(N,IEGEN,IORD) C C A=MATRIX TO BE DIAGONALIZED C N=DIMENSION OF SUBMATRIX TO BE DIAGONALIZED C IEGEN=0 EIGENVALUES AND EIGENVECTORS C =1 EIGENVALUES ONLY C IORD=0 NO ORDERING OF EIGENVALUES OR EIGENVECTORS (ORDER IN = ORDER OUT) C =1 ORDERING BY SIZE OF EIGENVALUES C EIVR=TRANSFORMATION MATRIX (MATRIX OF EIGENVECTORS) C C THIS ROUTINE USES A VARIABLE THRESHOLD JACOBI METHOD C use globdat c IMPLICIT REAL(8) (A-H,O-Z) c write(*,'('' starting '',i4,''x'',i4)')n,n IF(N.GT.0)GOTO 1 T(1,1)=1.0 RETURN 1 IF(IEGEN.GT.0)GOTO 102 DO 101 J=1,N DO 100 I=1,N T(I,J)=0.0d0 100 CONTINUE ! Fortran2018 T(J,J)=1.0d0 101 CONTINUE ! Fortran2018 C C FIND THE ABSOLUTELY LARGEST ELEMENT OF A C 102 ATOP=0.0D0 DO 1111 I=1,N ! Fortran2018 DO 111 J=1,N IF(ATOP.GE.DABS(H1(I,J)))GOTO 111 ATOP=DABS(H1(I,J)) 111 CONTINUE 1111 CONTINUE ! Fortran2018 c IF(ATOP)109,109,113 IF(ATOP.GT.0.d0)goto 113 ! Fortran2018 RETURN ! old label 109 C C CALCULATE THE STOPPING CRITERION -- DSTOP C 113 AVGF=(N*(N-1))*0.55D0 D=0.0D0 DO 1114 JJ=2,N ! Fortran2018 DO 114 II=2,JJ S=H1(II-1,JJ)/ATOP D=S*S+D 114 CONTINUE 1114 CONTINUE ! Fortran2018 DSTOP=(1.0D-06)*D C C CALCULATE THE THRESHOLD -- THRSH C THRSH=DSQRT(D/AVGF)*ATOP C C START A SWEEP C 115 IFLAG=0 DO 1130 JCOL=2,N ! Fortran2018 JCOL1=JCOL-1 DO 130 IROW=1,JCOL1 AIJ=H1(IROW,JCOL) C C COMPARE THE OFF DIAGONAL ELEMENT WITH THRSH C IF(DABS(AIJ).LE.THRSH)GOTO 130 AII=H1(IROW,IROW) AJJ=H1(JCOL,JCOL) S=AJJ-AII C C CHECK TO SEE IF THE CHOSEN ROTATION IS LESS THAN THE ROUNDING ERROR C IF SO, THEN DO NOT ROTATE C IF(DABS(AIJ).LE.(1.0D-09*DABS(S)))GOTO 130 IFLAG=1 C C IF THE ROTATION IS VERY CLOSE TO 45 DEGREES, SET SIN AND COS C TO 1/(ROOT 2) C IF((1.0D-10*DABS(AIJ)).LT.DABS(S))GOTO 116 S=0.70710678118654752D0 C=S GOTO 120 C C CALCULATE SIN AND COS FOR ROTATION THAT IS NOT VERY CLOSE C TO 45 DEGREES C 116 TT=AIJ/S S=0.25D0/DSQRT(0.25D0+TT*TT) C C COS=C, SIN=S C C=DSQRT(0.5D0+S) S=2.0D0*TT*S/C C C CALCULATION OF THE NEW ELEMENTS OF MATRIX A C 120 DO 121 I=1,IROW TT=H1(I,IROW) U=H1(I,JCOL) H1(I,IROW)=C*TT-S*U H1(I,JCOL)=S*TT+C*U 121 CONTINUE ! Fortran2018 I2=IROW+2 IF(I2.GT.JCOL)GOTO 123 CONTINUE DO 122 I=I2,JCOL TT=H1(I-1,JCOL) U=H1(IROW,I-1) H1(I-1,JCOL)=S*U+C*TT H1(IROW,I-1)=C*U-S*TT 122 CONTINUE ! Fortran2018 123 H1(JCOL,JCOL)=S*AIJ+C*AJJ H1(IROW,IROW)=C*H1(IROW,IROW)-S*(C*AIJ-S*AJJ) DO 124 J=JCOL,N TT=H1(IROW,J) U=H1(JCOL,J) H1(IROW,J)=C*TT-S*U H1(JCOL,J)=S*TT+C*U 124 CONTINUE C C ROTATION COMPLETED C SEE IF EIGENVECTORS REQUIRED C IF(IEGEN.GT.0)GOTO 126 DO 125 I=1,N TT=T(I,IROW) T(I,IROW)=SNGL(C*TT-T(I,JCOL)*S) T(I,JCOL)=SNGL(S*TT+T(I,JCOL)*C) 125 CONTINUE ! Fortran2018 C C CALCULATE THE NEW NORM D AND COMPARE WITH DSTOP C 126 CONTINUE S=AIJ/ATOP D=D-S*S IF(D.GE.DSTOP)GOTO 129 C C RECALCULATE DSTOP AND THRSH TO DISCARD ROUNDING ERRORS C D=0.0D0 DO 1128 JJ=2,N ! Fortran2018 DO 128 II=2,JJ S=H1(II-1,JJ)/ATOP D=S*S+D 128 CONTINUE 1128 CONTINUE ! Fortran2018 DSTOP=(1.0D-06)*D 129 THRSH=DSQRT(D/AVGF)*ATOP 130 CONTINUE 1130 CONTINUE ! Fortran2018 C IF(IFLAG.NE.0)GOTO 115 C C ARRANGE THE EIGENVALUES IN THE ORDER OF INCREASING ENERGY C ARRANGE THE EIGENVECTORS IN THE SAME ORDER C IF(IORD.EQ.0)RETURN NU=N DO 11 I=1,N IF(I.GE.NU)RETURN AMIN=H1(I,I) DO 10 J=I,NU IF(H1(J,J).GE.AMIN)GOTO 10 C C IF IEGEN IS -1 EXCLUDE UNCONVERGED EIGENVALUES FROM THIS ORDERING C TE=ABS(T(N,J))+ABS(T(N-1,J)) IF((TE.GT.0.05D0).AND.(IEGEN.EQ.-1))GOTO 15 II=I AMIN=H1(J,J) H1(J,J)=H1(I,I) H1(I,I)=AMIN 16 DO 12 K=1,N TEMP=T(K,II) T(K,II)=T(K,J) T(K,J)=SNGL(TEMP) 12 CONTINUE ! Fortran2018 GOTO 10 15 AM=H1(J,J) H1(J,J)=H1(NU,NU) H1(NU,NU)=AM II=NU NU=NU-1 GOTO 16 10 CONTINUE 11 CONTINUE c RETURN END c c_____________________________________________________________________________ C SUBROUTINE INV(N) C C INVERSION OF A MATRIX (Basically SSP vintage) C PARAMETER (NCONST=36,NSQCON=NCONST*NCONST) COMMON /BLKINV/AA REAL(8) AA(NCONST,NCONST),A(NSQCON) INTEGER(2) L(NCONST),M(NCONST) REAL(8) D,BIGA,HOLD C IJ=0 DO 6000 J=1,N DO 5000 I=1,N IJ=IJ+1 A(IJ)=AA(I,J) 5000 CONTINUE 6000 CONTINUE ! Fortran2018 C D=1.0D0 NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 2000 J=K,N ! Fortran2018 IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I c10 IF(DABS(BIGA)-DABS(A(IJ))) 15,20,20 IF(DABS(BIGA)-DABS(A(IJ)).ge.0.d0)goto 20 ! Fortran2018 old label 10 BIGA=A(IJ) ! old label 15 L(K)=I M(K)=J 20 CONTINUE 2000 CONTINUE ! Fortran2018 C J=L(K) c IF(J-K) 35,35,25 IF(J-K.LE.0)goto 35 ! Fortran2018 KI=K-N ! old label 25 DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) A(JI)=HOLD 30 CONTINUE C 35 I=M(K) C IF(I-K) 45,45,38 IF(I-K.LE.0)goto 45 ! Fortran2018 JP=N*(I-1) ! old label 38 DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) A(JI)=HOLD 40 CONTINUE C c45 IF(BIGA) 48,46,48 45 IF(BIGA.NE.0.d0)goto 48 ! Fortran2018 D=0.0D0 ! old label 46 WRITE(*,'(1X/'' ***** ZERO DETERMINANT *****''/)') RETURN 48 DO 55 I=1,N c IF(I-K) 50,55,50 IF(I-K.EQ.0)goto 55 ! Fortran2018 IK=NK+I ! old label 50 A(IK)=A(IK)/(-BIGA) 55 CONTINUE C DO 1065 I=1,N ! Fortran2018 IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N c IF(I-K) 60,65,60 IF(I-K.EQ.0)goto 65 ! Fortran2018 c60 IF(J-K) 62,65,62 IF(J-K.EQ.0)goto 65 ! Fortran2018 old label 60 KJ=IJ-I+K ! old label 62 A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE 1065 CONTINUE ! Fortran2018 C KJ=K-N DO 75 J=1,N KJ=KJ+N C IF(J-K)70,75,70 IF(J-K.EQ.0)goto 75 ! Fortran2018 A(KJ)=A(KJ)/BIGA ! old label 70 75 CONTINUE C D=D*BIGA C A(KK)=1./BIGA 80 CONTINUE C K=N 100 K=K-1 c IF(K) 150,150,105 IF(K.LE.0)goto 150 ! Fortran2018 I=L(K) ! old label 105 c IF(I-K) 120,120,108 IF(I-K.LE.0)goto 120 ! Fortran2018 JQ=N*(K-1) ! old label 108 JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) A(JI)=HOLD 110 CONTINUE ! Fortran2018 120 J=M(K) c IF(J-K) 100,100,125 IF(J-K.le.0)goto 100 ! Fortran2018 KI=K-N ! old label 125 DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) A(JI)=HOLD 130 CONTINUE GOTO 100 C 150 DO 6001 I=1,N DO 5001 J=1,N IJ=N*(J-1)+I AA(I,J)=A(IJ) 5001 CONTINUE 6001 CONTINUE C RETURN END c c_____________________________________________________________________________ C SUBROUTINE SORTC(N,M) IMPLICIT INTEGER (I-N) PARAMETER (NLIN=2000,msort=10000) COMMON /sortcc/wk,IPT INTEGER(4) IPT(msort) REAL(8) WK(msort),EE C C ... This routine sorts the quantities in vector WK from N to M in C ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities C if(M.le.N) return DO 101 I=N,M-1 J=I 106 J=J+1 c IF(WK(J)-WK(I))103,104,104 IF(WK(J)-WK(I).GE.0)goto 104 ! Fortran2018 EE=WK(I) ! old label 103 WK(I)=WK(J) WK(J)=EE K=IPT(I) IPT(I)=IPT(J) IPT(J)=K 104 IF(J.EQ.M)GOTO 101 GOTO 106 101 CONTINUE C RETURN END C c_____________________________________________________________________________ C SUBROUTINE SORTJ(JJ,N,M) IMPLICIT INTEGER (I-N) PARAMETER (NLIN=2000,msort=10000) COMMON /sortcc/wk,IPT INTEGER(4) IPT(msort) REAL(8) WK(msort),eout(msort),ee C C ... This routine sorts the quantities in vector WK from N to M which are C assumed to correspond to the complete set of eigenvalues for a J-block c defined by J C It is assumed that eigenvalues are in order K=...-2,-1,0,1,2... and C are sorted according to ascending value of Tau C j=jj/2 do 1 k=0,J if(K.eq.0)then indx=0 eout(n+indx)=wk(n+j) ipt(n+indx)=n+j else indx=indx+1 eout(n+indx)=wk(n+j-k) ipt(n+indx)=n+j-k indx=indx+1 eout(n+indx)=wk(n+j+k) ipt(n+indx)=n+j+k if(eout(n+indx).lt.eout(n+indx-1))then ee=eout(n+indx-1) eout(n+indx-1)=eout(n+indx) eout(n+indx)=ee ii=ipt(n+indx-1) ipt(n+indx-1)=ipt(n+indx) ipt(n+indx)=ii endif endif 1 continue c do 2 k=n,m wk(k)=eout(k) 2 continue C RETURN END C c_____________________________________________________________________________ C C Actual length of a string (equivalent to LEN_TRIM) C function lentrm(carg) character carg*(*) c nn=len(carg) do 1 n=nn,1,-1 if(carg(n:n).gt.' ')goto 2 1 continue 2 lentrm=n c return end c_____________________________________________________________________________ C_____________________________________________________________________________ C C Modules used by QS: C C globdat - specification of variables and arrays for global use C gleout - gle output C confor - constant and error formatting for output C read - to set up the Hamiltonian C mhigh - establish dimensioning for H C quantn - return quantum numbers for an index in H C numenl - the index in H for given quantum numbers C setham - used by READ to set up H matrix c hamil - matrix element in H C threej - evaluate the 3J symbol C sixj - evaluate the 6J symbol C x9j - evaluate the 9J symbol C triang - check the triangular condition C tri C fac C suml C sumk C mint - integer/half-integer check used by X9J C lsqrw - weighted least squares C hdiagm - jacobi matrix diagonalisation c inv - matrix inversion C sortc - magnitude sort c sortj - sort according to increasing tau C lentrm - substitute for LEN_TRIM C C_____________________________________________________________________________ C_____________________________________________________________________________