C---------------------------------------------------------------------- C Fitting of Magnetic susceptibility as a function of temperature. C Ref: G. V. R. Chandramouli, C. Balagopalakrishna C M. V. Rajasekharan, and P. T. Manoharan C Computers Chem. 20 (1996), 353-358. C Date: September, 1994, modified April,'02. C C ---------------------------------------------------------------------- C C The maximum number of data points are 300 C The maximum number of minimisation parameters are 12 C The maximum number of terms in function are 10 C The input data file contains C NP - number of points C T(I) - Temperature, XM(I) - Observed Susceptibility per ion C at T(I) including diamagnetic correction. I runs from 1 to NP. C The remaining data is given on terminal. C g-value can be fixed or treated as a minimization parameter. C Functions considered are: C 1. Curie 4. Bleaney-Bower 7. Modified BF 10. DIM C 2. Zero-field 5. Magnetization 8. Dimers C 3. Ising 6. Bonner-Fisher 9. Other C Other: 1. User defined 3. Hatfield's tetramer (S=1/2) C 2. Trimer (S=1/2) 4. Kahn's dimer chains C The user defined function requires modification of subroutines C USRPAR and USER. See comments in these routines. C Interdimer interactions are included by molecular field approx. C Weiss term is optionally included for Curie term. C One can combine the above functions (terms) to choose an expression C for the susceptibility calculation. A menu of preselected C combinations is presented to opt the desired expression. C Also there is a provision for other combinations (option 99). C C Initially a plot of observed and calculated data is shown. C User can choose one of the three plot options -- C SPLOT for graphics screen. PPLOT -- for text screen C HPPLOT to create an output file with HP plot instructions C The initial parameters of the expression may be interactively C modified. C Output data is written on SUSCEP.OUT only if iterations were done. C SUSCEP.HP file is created only HPPLOT option is chosen. C EXPLMT is machine dependent limit of exponent of base e (=2.714) C C Convention of J value: The singlet-triplet gap is treated as 2J. C Unit of J is wave number. Calculates susceptibility per magnetic ion. C C --------------------------------------------------------------------- C Program uses Microsoft 5.0 graphics routines for screen graphics C Tested with Microsoft FORTRAN compiler ver. 5.0 C The subroutine NELMIN is picked from STATISTICAL ALGORITHMS 47. C Minor modifications were made to this to extract error estimates. C Error Estimation is done according to C Gregory R. Phillips and Edward M. Eyring C Anal. Chem. 1988, vol 60, pp 738-741 and references therein. C The routine GJRD was taken from SCIENTIFIC SUBROUTINE PACKAGE C --------------------------------------------------------------------- INCLUDE 'FGRAPH.FI' C --------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) CHARACTER PAR*8, INFILE*20, OPTION*1 DOUBLE PRECISION MIN(12), TC(300), XIC(300),XMO(300) DIMENSION SVPAR(12),SVPEP(12),SIGMA(12) COMMON /BLK1/ T(300), XM(300), DA, NP, IENTRY COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK2A/ BK, BN, GBK, GBN, TIP COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK4/ START(12), STEP(12), PAR(12) COMMON /BLK5/ IN, IO COMMON /BLK6/ ITYPE, NFN, NF(10), NPAR COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT C Input section (also look subroutine GETPAR) IN = 5 IO = 6 EXPLMT = DLOG(1.0D307) NPRLMT=12 NFNLMT=10 C The observed data points are read from INFILE WRITE(IO,5) 5 FORMAT(1H1,///T27,' SUSCEPTIBILITY DATA FITTING '/ 1 T22,' by G.V.R.C, C.B.G., M.V.R., & P.T.M'/ 2 T22,' Ref: Computers & Chemistry '/// 3 ' Experimental data file name?') READ(IN,10) INFILE 10 FORMAT(A) OPEN(UNIT=1,FILE=INFILE,STATUS='OLD') READ(1,*) NP READ(1,*) (T(I), XMO(I), I=1,NP) WRITE(IO,*)'No. of observed data points:', NP WRITE(IO,*) 1 'Temp.(K) X(M) Temp.(K) X(M) Temp.(K) X(M)' WRITE(IO,11)(T(I),XMO(I), I=1,NP) 11 FORMAT(3(F7.2,F9.6,5X)) CLOSE(1) WRITE(IO,*)'If the above X(M) values are NOT per magnetic ion,' WRITE(IO,*)'multiply them by a Factor to get X(M) value per ion' WRITE(IO,*)'See next menu. Press Enter key...' READ(IN,10) OPTION C Default settings DIA=0.D0 TIP0=0.D0 FACTOR=1.D0 NTH=0 NEQN=24 G=0.D0 12 NPAR = 0 C Present Main Menu WRITE(IO,20) DIA,TIP0 20 FORMAT(' A. Diamagnetic value to be subtracted from observed' 2 ' data=',1PE10.3/' B. TIP value to be added to the calculated', 3 ' Susceptibility=',1PE10.3) WRITE(IO,21) FACTOR,NEQN 21 FORMAT(' C. Multipy the Observed data by factor=',F5.2,/ 5 ' D. Show the present values'/ 6 ' E. Expression No. =',I2,/ 7 ' (A two digit no. MN (01-99). First digit M selects the ', 8 'type of expression.'/ 9 ' Second digit N selects the function F(N). ', 9 'MN=99 for a new expression)') WRITE(IO,25) 25 FORMAT( 1X,75('-')/' M Expressions for X(cal)' 1 ,T35,'| N F(N)',T55,'Consts', T62,'Variables'/ 2 1X,75('-')/ 3 ' 0 Simple F(N)',T35,'| 1. Curie', T55, 'S'/ 4 ' 1 F(N) + DIM',T35,'| 2. Zero-Field', T55, 'S', T62, 'D'/ 5 ' 2 x F(N) + (1-x) Curie',T35,'| 3. Ising', T62, 'J-ISNG'/ 6 ' 3 x F(N) + (1-x) Z-F',T35,'| 4. Bleaney-Bower',T62,'J-BB'/ 7 ' 4 x F(N) + (1-x) Curie + DIM',T35, 8 '| 5. Magnetization',T55,'B',T62,'J-MAGN') WRITE(IO,30) G 30 FORMAT(' 5 x F(N) + (1-x) Z-F + DIM',T35, 1 '| 6. Bonner-Fisher', T62, 'J-BF'/ 2 ' 9 None of the above',T35, 3 '| 7. Modified BF', T62, 'J-MBF, Alpha'/ 4 ' x is COEFF1 (<1.0)',T35, 5 '| 8. Dimers', T55, 'S1, S2', T62, 'J-Dimer'/ 6 ' e.g. No.24 is Bleaney-Bower with',T35,'| 9. Other functions'/ 7 ' Curie impurity',T35,'|10. Constant (DIM)', T62, 'DIM'/ 8 1X,75('-')/' F. Proceed for calculation'/ 9 ' G. g-value (g=0 to make it fitting parameter)=',F7.4) 31 WRITE(IO,*)'To change any value, ENTER OPTION (A-G):' READ(IN,10) OPTION IF(OPTION.EQ.'A' .OR.OPTION.EQ.'a') GOTO 32 IF(OPTION.EQ.'B' .OR.OPTION.EQ.'b') GOTO 34 IF(OPTION.EQ.'C' .OR.OPTION.EQ.'c') GOTO 36 IF(OPTION.EQ.'D' .OR.OPTION.EQ.'d') GOTO 38 IF(OPTION.EQ.'E' .OR.OPTION.EQ.'e') GOTO 40 IF(OPTION.EQ.'F' .OR.OPTION.EQ.'f') GOTO 46 IF(OPTION.EQ.'G' .OR.OPTION.EQ.'g') GOTO 42 GOTO 31 32 WRITE(IO,*)'Enter diamagnetic correction including -ve sign' READ(IN,*) DIA GOTO 31 34 WRITE(IO,*)'Enter Temperature Independent Paramagnetism value' READ(IN,*) TIP0 GOTO 31 36 WRITE(IO,*)'Enter the factor' READ(IN,*) FACTOR GOTO 31 38 GOTO 12 40 WRITE(IO,*)'Enter S.No. (1-99) for equation' READ(IN,*) NEQN GOTO 31 42 WRITE(IO,*)'Enter g-value (g=0 to make it fit parameter)' READ(IN,*) G GOTO 31 C Proceed for calculation 46 CONTINUE DO 37 I=1,NP 37 XM(I)=(XMO(I)-DIA)*FACTOR TIP=TIP0*FACTOR DA=0.D0 DO 48 I=1,NP 48 DA=DA+XM(I)*XM(I) C The g-value below 1.E-6 makes g as minimization parameter. IF(DABS(G).LT.1.D-6) CALL GETPAR('g-value ') IF (NEQN/10.EQ.9) GOTO 200 C Separate the two digits from NEQN to select the expression J = NEQN - NEQN/10*10 C 0, 10, 20, 30, 40, 50, 60, 70, 80-98 are unused. Skip them. C NFN is number of terms in expression. See later for ITYPE IF (J.EQ.0 .OR. NEQN .GT. 79) GOTO 40 I = NEQN/10 ITYPE = 1 NFN = 1 NF(1) = J IF (I.GT.0) GOTO 50 GOTO 300 50 NFN = 2 IF (I .GT. 3) GOTO 100 IF (I .GT. 1) GOTO 60 ITYPE = 1 NF(2) = 10 GOTO 300 60 ITYPE = 2 GOTO (50,70,80), I 70 NF(2) = 1 GOTO 300 80 NF(2) = 2 GOTO 300 100 NFN = 3 NF(3) = 10 I = I -2 GOTO 60 C Expression is defined by user by giving NFN, ITYPE C ITYPE =1, 2, or 3 as described in FORMAT label 220 C A, B, C,... are the Functions and x = COEFF1, y = COEFF2 C 200 WRITE(IO,205) 205 FORMAT(' List of Functions:'/' 1. Curie Law'/ 1 ' 2. Zero-Field'/' 3. Ising'/' 4. Bleaney-Bower'/ 2 ' 5. Magnetization'/' 6. Bonner-Fisher'/ 3 ' 7. Modified BF'/' 8. Dimers'/' 9. Other functions'/ 4 ' 10. Constant (DIM)') WRITE(IO,*)'Enter No. of Terms (Fns.) in Eqn (Max. ',NFNLMT,')' READ(IN,*) NFN WRITE(IO,*)'Enter S. Nos. of Fns. from above list' READ(IN,*) (NF(I), I=1,NFN) 210 WRITE(IO,220) 220 FORMAT(' ITYPE Type of Equation'/ 1 ' 1. A + B + C + ...'/ 2 ' 2. x(A) + (1-x)(B) + C + ...'/ 3 ' 3. y [ x(A) + (1-x)(B) ] + (1-y) (C) + D + ...'/ 4 ' Enter ITYPE (1-3)') READ (IN,*) ITYPE IF (ITYPE .LT. 1 .OR. ITYPE .GT. 3) GOTO 210 C Define the initial guess values for minimization parameters 300 DO 310 I = 1,NFN CALL DEFPAR( NF(I), I) 310 CONTINUE IF(ITYPE.EQ.1) GOTO 400 CALL GETPAR ('COEFF1 ') IF(ITYPE.EQ.3) CALL GETPAR('COEFF2 ') C Proceed for calculation of X(cald) C constants used through out the program C HCBIK = 6.62608 * 2.9979246 / 1.38066 / 10 400 HCBIK = 1.4387676 C B2NBIK = 9.27402**2 * 6.02214 / 1.38066 / 1000 B2NBIK = 0.37514586 C HCNB2 = hc/Nb2 6.62608*2.9979246/(6.02214*9.27402**2)*100 HCNB2 = 3.8352217 C BK = beta*k BN=beta*N BK = 6.7170918D-5 BN = 5584.9447D0 IF (G.LE.1.D-6) GOTO 405 HCNGB2 = HCNB2/G/G G2B2NK = G*G*B2NBIK GBK = G*BK GBN = G*BN C Select type of plot and prepare plot data 405 WRITE(IO,410) 410 FORMAT(/' Plot options:'/' 1. X vs. T'/' 2. Mu vs. T'/ 1 ' 3. XT vs. T'/' 4. 1/X vs. T'/' 5. X vs 1/T'/' 6. QUIT' 2 /' Enter Option (1-6) ') READ(IN,*) NPLOP WRITE(IO,*) 'Plot type: 1. Splot/ 2.Pplot/ 3.HPPlot' WRITE(IO,*) 'Enter Option (1-3)' READ(IN,*) IPLOP CALL XDATA(T, NP,TC,NC,NPLOP,IPLOP) 415 DO 420 I=1,NPAR 420 MIN(I) = START(I) DO 430 I=1,NC 430 XIC(I) = FNO(MIN,TC(I)) CALL YDATA(T,XM,NP,TC,XIC,NC,NPLOP) IF(IPLOP.EQ.1) CALL SPLOT(NP,NC) IF(IPLOP.EQ.2) CALL PPLOT(NP) IF(IPLOP.EQ.3) CALL HPPLOT(NP,NC,NPLOP) C Interactively change the initial values of parametes WRITE(IO,*) 'No. Parameter Value' WRITE(IO,440) (I,PAR(I),MIN(I),I=1,NPAR) 440 FORMAT(I3,'. ',A,F12.6) 450 WRITE(IO,*) 1 'Change Parameter No. (0-Plot, -1 Iteration, -2 End):' READ(IN,*) I IF (I.GT.NPAR .OR. I.LT.-2) GOTO 450 IF (I.EQ.0) GOTO 415 IF (I.EQ.-1) GOTO 460 IF (I.EQ.-2) GOTO 700 WRITE(IO,*) 'Enter New value of ',PAR(I) READ(IN,*) START(I) GOTO 450 C Happy with initial values of minimization parameters C go ahead with iteration 460 REQMIN = 1.0 D-16 KONGVE = 1 IENTRY = 0 ICOUNT=200 DO 470 I=1,NPAR SIGMA(I)=0.0D0 STEP(I)=START(I)*1.D-2 SVPEP(I)=STEP(I) 470 SVPAR(I)=START(I) WRITE(IO,490) (PAR(I),I=1,NPAR),'L.S.Error' YNEWLO=FN(START) 480 WRITE(IO,*)'No. of iterations',ICOUNT WRITE(IO,481)REQMIN 481 FORMAT(' Convergence limit for iteration:',1PE10.2) WRITE(IO,*)' Parameter Value Step' DO 482 I=1,NPAR 482 WRITE(IO,483) I,PAR(I), START(I), STEP(I) 483 FORMAT(1X,I2,'. ',A,2F12.6) IF(ICOUNT.EQ.0) GOTO 500 WRITE(IO,*) 1 'Change step No. (0 to proceed, -1 No. of Iters. ', 2 '-2 Convergence):' READ(IN,*) I IF (I.EQ.-2) GOTO 485 IF (I.EQ.-1) GOTO 487 IF (I.EQ.0) GOTO 488 WRITE(IO,*) 'Enter New step value of ',PAR(I) READ(IN,*) STEP(I) GOTO 480 485 WRITE(IO,*) 'Enter REQMIN value for convergence:' READ(IN,*) REQMIN GOTO 480 487 WRITE(IO,*)'Enter No. of iterations' READ(IN,*) ICOUNT GOTO 480 488 DO 489 I=1,NPAR 489 SVPEP(I)=STEP(I) WRITE(IO,490) (PAR(I),I=1,NPAR),'L.S.Error' 490 FORMAT(8X,7(1X,A,1X)) CALL NELMIN (NPAR, START, MIN, YNEWLO, REQMIN, STEP, 1 KONGVE, ICOUNT) WRITE(*,*)'Iteration is completed. Calculation ESDs' CALL ERREST(MIN, SIGMA,NPAR,NP,IER) WRITE(IO,490) (PAR(I),I=1,NPAR),'L.S.Error' YNEWLO=FN(MIN) C Output the results, also write on file XIFIT.OUT 500 OPEN (UNIT=2, FILE='SUSCEP.OUT', STATUS='UNKNOWN') CALL EQN WRITE(IO,510) WRITE(2,510) 510 FORMAT(//' Starting:'/' Parameter Value Step') WRITE(IO,520) (PAR(I),SVPAR(I), SVPEP(I), I=1,NPAR) WRITE(2,520) (PAR(I),SVPAR(I), SVPEP(I), I=1,NPAR) 520 FORMAT(1X, A, 1PE12.5,5X,1PE12.5) WRITE(IO,540) (PAR(I), MIN(I),SIGMA(I), I=1,NPAR) WRITE(2,540) (PAR(I), MIN(I),SIGMA(I), I=1,NPAR) 540 FORMAT(/' Solution:'/' Parameter Value Std. devn.'/ 1 (1X, A,' =',1PE12.5,' +/-',1PE12.4)) IF(IER.EQ.1) WRITE(IO,545) IF(IER.EQ.1) WRITE(2,545) 545 FORMAT(' Sigma values could not be estimated correctly.') IF (DABS(G) .LT.1.E-6) GOTO 555 WRITE(IO,550) G WRITE(2,550) G 550 FORMAT(' g-value given =', F8.4) 555 WRITE(IO,560) YNEWLO WRITE(2,560) YNEWLO 560 FORMAT(' L.S. Error =', D12.4/) C Generate L. S. Fitted values WRITE(IO,*) 'Press RETURN key' READ(IN,*) WRITE(IO,570) WRITE(2,570) 570 FORMAT(4X, 'Temp.(K)', 8X, 'X-Expt', 9X, 'X-Cald.', 1 9X,'Difference') DO 590 I=1,NP XIC(I) = FNO(MIN,T(I)) DIFF = XM(I)-XIC(I) WRITE(IO,580) T(I), XM(I), XIC(I), DIFF WRITE(2,580) T(I), XM(I), XIC(I), DIFF 580 FORMAT(F15.2,3F15.5) 590 CONTINUE WRITE(2,*)'Diamagnetic correction made on X-obsd =',DIA WRITE(2,*)'TIP added to X-cald. =',TIP WRITE(2,*)' TEMP. X-CALD.' XINT=(T(NP)-T(1))/300 DO 600 I=1,301 TEMP=(I-1)*XINT+T(1) XI = FNO(MIN,TEMP) 600 WRITE(2,610) TEMP, XI 610 FORMAT(1X,F7.2,1X,1PE13.6) CLOSE(2) C The iteration results are saved for future use DO 620 I=1,NPAR 620 START(I) = MIN(I) GOTO 415 700 WRITE(IO,710) 710 FORMAT(' Options:'/' 1. Select Function'/' 2. Plot option'/ 1 ' 3. End'/ ' Enter option (1-3)') READ(IN,*) NREPT GOTO (12,405,999), NREPT 999 STOP END C Defvar routine defines the parameters which are to be minimised. SUBROUTINE DEFPAR(NF,I) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER ANS*1 COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK5/ IN, IO COMMON /BLK9/ NOTH COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT IF(I.LT.NFNLMT .AND. I.GT.0) GOTO 1 WRITE(IO,*) 'No. of terms cannot be more than ',NFNLMT STOP 1 IF(NF.GT.0 .AND. NF.LE.10) GOTO 2 WRITE(IO,*)'Function not available in the list (1-10)' STOP 2 IF(NF.GT.2) GOTO 30 GOTO (10,20), NF 10 WRITE(IO,*) 'Enter Spin (S) for Curie term' READ(IN,*) S(I) GOTO 150 20 WRITE(IO,*)'Zero-field Equations:' WRITE(IO,*)'Enter Spin (1.0, 1.5, 2.0, or 2.5)' READ(IN,*) S(I) NSI=S(I)*2-0.9 IF( NSI.LT.1 .OR. NSI.GT.4) GOTO 30 IF( NSI.EQ.1) CALL GETPAR('D-1 ') IF( NSI.EQ.2) CALL GETPAR('D-3/2 ') IF( NSI.EQ.3) CALL GETPAR('D-2 ') IF( NSI.EQ.4) CALL GETPAR('D-5/2 ') GOTO 150 30 NFM2=NF-2 GOTO (40,50,60,70,80,90,110,100), NFM2 40 WRITE(IO,*) 'Ising Eqn. parameters' CALL GETPAR('J-Ising ') GOTO 150 50 WRITE(IO,*) 'Bleaney-Bower Eqn. parameters' CALL GETPAR('J-BB ') GOTO 150 60 WRITE(IO,*) 'Magnetization eqn. parameters' CALL GETPAR('J-MAGN ') WRITE(IO,*)'Enter Field (B)' READ(IN,*) Z(I) GOTO 150 70 WRITE(IO,*) 'Bonner Fisher Eqn. parameters' CALL GETPAR('J-BF ') GOTO 150 80 WRITE(IO,*) 'Modified Bonner Fisher Eqn. parameters' CALL GETPAR('J-MBF ') CALL GETPAR('ALFA-MBF') GOTO 150 90 WRITE(IO,*) 'Dimer Eqn. parameters' WRITE(IO,*)'Enter S1, S2 of dimer' READ(IN,*) S(I),Z(I) CALL GETPAR('J-Dimer ') GOTO 150 100 CALL GETPAR('DIM ') NS(I)=0 GOTO 160 110 CALL OTHPAR(NF,I) 150 IF(NF.GT.1) 1 WRITE(IO,*)'Include Neighboring interaction Z*JP (Y/N)?' IF(NF.EQ.1) WRITE(IO,*)'Include Weiss constant (theta) (Y/N)?' READ(IN,'(A)') ANS NS(I)=0 IF (.NOT.(ANS.EQ.'Y' .OR. ANS.EQ.'y')) GOTO 160 NS(I)=1 IF(NF.EQ.1)CALL GETPAR('THETA ') IF(NF.GT.1)CALL GETPAR('Z*JP ') 160 RETURN END C GETPAR reads the values of parameters to be minimised and assigns C number of variables NPAR SUBROUTINE GETPAR(PARNAM) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*8 PARNAM, PAR COMMON /BLK4/ START(12), STEP(12), PAR(12) COMMON /BLK5/ IN, IO COMMON /BLK6/ ITYPE, NFN, NF(10), NPAR COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT NPAR = NPAR+1 IF(NPAR .LE.NPRLMT) GOTO 10 WRITE(IO,*) 'Minimisation parameters should not exceed ',NPRLMT STOP 10 PAR(NPAR)=PARNAM WRITE(IO,*) 'Enter Initial value for ',PAR(NPAR) READ(IN,*) START(NPAR) RETURN END C Calculation of Total error - Supplies to NELMIN DOUBLE PRECISION FUNCTION FN(A) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK1/ T(300), XM(300), DA, NP, IENTRY COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK5/ IN, IO COMMON /BLK6/ ITYPE, NFN, NF(10), NPAR DIMENSION A(1) IENTRY=IENTRY+1 DIF = 0.D0 DO 10 I=1,NP FC=FNO(A,T(I)) X = XM(I)-FC DIF = DIF +X*X 10 CONTINUE FN = DIF/DA WRITE(IO,20) IENTRY, (A(I),I=1,NPAR), FN 20 FORMAT(I5,(7F10.4)) RETURN END C FNO evaluates the value of expression. The calculated C susceptibility of dimer, trimer, etc. are divided by 2, 3, etc. C to get susceptibility per magnetic ion. DOUBLE PRECISION FUNCTION FNO(A,T) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK2A/ BK, BN, GBK, GBN, TIP COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK5/ IN, IO COMMON /BLK6/ ITYPE, NFN, NF(10), NPAR COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT DIMENSION A(1), F(10) IPAR = 0 IF (DABS(G) .GT. 1.D-6) GOTO 1 IPAR = IPAR + 1 G2B2NK = A(1)*A(1)*B2NBIK HCNGB2 = HCNB2/A(1)/A(1) GBK = A(1)*BK GBN = A(1)*BN 1 IF(NFN.LE.NFNLMT) GOTO 5 WRITE(IO,*) 'Number of terms exceed ',NFNLMT STOP 5 DO 105 I=1,NFN N = NF(I) GOTO (10, 20, 30, 40, 50, 60, 70, 80, 95, 90), N 10 CALL CURIE(S(I), T, F(I)) GO TO 100 20 IPAR = IPAR+1 CALL ZF(S(I), A(IPAR), T, F(I)) GOTO 100 30 IPAR = IPAR+1 CALL ISING(A(IPAR), T, F(I)) F(I)=F(I)/2.D0 GOTO 100 40 IPAR = IPAR+1 CALL BB(A(IPAR), T, F(I)) F(I)=F(I)/2.D0 GOTO 100 50 IPAR = IPAR+1 GG=G IF (DABS(G) .GT. 1.D-6) GG=A(1) CALL MAGN(A(IPAR), Z(I), T, F(I)) GOTO 100 60 IPAR = IPAR+1 CALL BF(A(IPAR), T, F(I)) GOTO 100 70 IPAR = IPAR+2 CALL MBF(A(IPAR-1), A(IPAR), T, F(I)) F(I)=F(I) GOTO 100 80 IPAR = IPAR+1 CALL DIMER (S(I), Z(I), A(IPAR), T, F(I)) F(I)=F(I)/2.D0 GOTO 100 90 IPAR = IPAR+1 F(I) = A(IPAR) GOTO 100 95 CALL OTHER(A, IPAR, T, F(I),I) 100 IF(.NOT.(NS(I).EQ.1)) GOTO 105 IPAR = IPAR+1 IF(N.EQ.1.AND.A(IPAR).NE.T) F(I)=F(I)*T/(T-A(IPAR)) C Note that this expression is also modified to consider C monomeric value. IF(N.GT.1) F(I) = F(I)/(1.D0-4.D0*HCNGB2*A(IPAR)*F(I)) 105 CONTINUE IF(ITYPE .NE.2) GOTO 110 F(1)=F(1)*A(NPAR) F(2)=F(2)*(1.D0-A(NPAR)) 110 IF(ITYPE.NE.3) GOTO 120 F(1) = F(1)*A(NPAR)*A(NPAR-1) F(2) = F(2)*A(NPAR)*(1.D0-A(NPAR-1)) F(3) = F(3)*(1.D0 - A(NPAR)) 120 FNO = 0.D0 DO 130 I=1,NFN 130 FNO = FNO+F(I) FNO=FNO+TIP RETURN END C Curie function SUBROUTINE CURIE(S,T,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO F = G2B2NK/T*S*(S+1.D0)/3.D0 RETURN END C S=1,3/2,2,5/2 functions with Zero-Field splitting SUBROUTINE ZF(S,D,T,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT NS=S*2-0.9 IF(NS.LT.1 .OR. NS.GT.4) GOTO 200 CON = G2B2NK/T/3.D0 IF (D.EQ.0.D0) GOTO 100 X=D*HCBIK/T GOTO (10,20,30,40), NS 10 IF(DABS(X).GT.EXPLMT) WRITE(IO,60) X EX = DEXP(-X) TOP1 =2.*EX TOP2 = (2.D0/X)*(1.D0+EX) BOT = 1.D0+2.D0*EX GOTO 50 20 X2=2.D0*X IF(DABS(X2).GT.EXPLMT) WRITE(IO,60) X2 EX2 = DEXP(-X2) TOP1=1.D0+9.D0*EX2 TOP2 = 4.D0 + (3.D0/X)*(1.D0-EX2) BOT = 4.D0*(1.D0+EX2) GOTO 50 30 X4 =4.D0*X IF(DABS(X4).GT.EXPLMT) WRITE(IO,60) X4 EX = DEXP(-X) EX4 = DEXP(-X4) TOP1 = 2.D0*EX+8.D0*EX4 TOP2 = (6.D0/X) *(1-EX)+(4.D0/3.0D0/X)*(EX-EX4) BOT = 1.D0+2.D0*(EX+EX4) GOTO 50 40 X6 = 6.D0*X IF(DABS(X6) .GT. EXPLMT) WRITE(IO,60) X6 EX2 = DEXP(-2.D0*X) EX6 = DEXP(-6.D0*X) TOP1 = 1.D0+9.D0*EX2 +25.D0*EX6 TOP2 = 9.D0+(8.D0/X)*(1.D0-EX2)+4.5D0/X*(EX2-EX6) BOT = 4.D0*(1.D0 + EX2 + EX6) 50 F = CON*(TOP1+TOP2+TOP2)/BOT GOTO 200 60 FORMAT(' Exponent Warning in Z-F routine-->',E12.5) 100 F = CON*S*(S+1) 200 RETURN END C Ising Eq. SUBROUTINE ISING(A,T,F) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT AJ=A*HCBIK/T IF (AJ.GT.EXPLMT) WRITE(IO,*) 1 'Exponent Warning in ISING: AJ = ',AJ X=DEXP(AJ) F=G2B2NK/(2.D0*T)*X/(1+X) RETURN END C Bleaney Bower's expression SUBROUTINE BB(AJ,T,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT X = 2.D0*HCBIK*AJ/T IF(DABS(X).GT.EXPLMT) + WRITE(IO,*)'Exponent Warning in BB--> X =',X TOP = 2.0*G2B2NK/T BOT = 3.D0+DEXP(-X) F = TOP/BOT RETURN END C Generalized magnetization eq. SUBROUTINE MAGN(AJ,B,T,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK2A/ BK, BN, GBK, GBN, TIP COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT X = -2.D0*HCBIK*AJ/T IF(DABS(X).GT.EXPLMT) + WRITE(IO,*)'Exponent Warning in MAGN--> X =',X TERM=GBK*B/T TOP=GBN*SINH(TERM) BOT=DEXP(X)+2.D0*COSH(TERM)+1.D0 F = TOP/BOT/B RETURN END C Bonner-fischer eq. SUBROUTINE BF(AJ,T,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO X=HCBIK*DABS(AJ)/T X2=X*X X3=X2*X ANUM = 0.25 + 0.149445*X + 0.30094*X2 DENO = 1 + 1.9862*X + 0.68854*X2 + 6.0626*X3 F = G2B2NK/T * ANUM/DENO RETURN END c Modified Bonner-Fisher equation c (this routine modified April,'02 -mvr) c Ref: W. E. Hatfield, J. Appl. Phys. 52 (1981),1985 and c "Molecular Magnetism", O. Kahn ..., p.263-265. SUBROUTINE MBF(AJ,AL,T,F) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO AL2=AL*AL AL3=AL2*AL AL4=AL3*AL IF(AL.LE.0.4) THEN A1=.25 B=-.12587+.22752*AL C=.019111-.13307*AL+.509*AL2-1.3167*AL3+1.0081*AL4 D=.10772+1.4192*AL E=-.0028521-.42346*AL+2.1953*AL2-.82412*AL3 F=.37754-6.702201E-2*AL+6.9805*AL2-21.678*AL3+15.838*AL4 ELSE A1=0.25 B=-.13695+.26387*AL C=.017025-.12668*AL+.49113*AL2-1.1977*AL3+.87257*AL4 D=.070509+1.3042*AL E=-.0035767-.40837*AL+3.4862*AL2-.73888*AL3 F=.36184-.065528*AL+6.65875*AL2-20.945*AL3+15.42504*AL4 ENDIF X=HCBIK*(-AJ)/T X2=X*X X3=X2*X ANUM = A1 + B*X +C*X2 DENO = 1 + D*X + E*X2 + F*X3 F = G2B2NK/T * ANUM/DENO RETURN END C Heisenberg-intradimer interaction for various spin systems SUBROUTINE DIMER (S1,S2,AJ,T,F) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT DIMENSION I(10), E(10) AJBIKT = HCBIK*AJ/T I1 = 2*S1+0.01 I2 = 2*S2+0.01 ITOT = I1+I2+1 IMIN = IABS(I1-I2)+1 K = 0 DO 40 LL = IMIN, ITOT, 2 L = LL-1 K = K+1 I(K) = L IF(LL.EQ.IMIN) IEK = 0 IF(LL.GT.IMIN) IEK = IEK + L TERM = IEK*AJBIKT IF(DABS(TERM).GT. EXPLMT) 1 WRITE(IO,*)'Exponent Warning in DIMER-->, TERM =',TERM E(K) = DEXP(TERM)*LL 40 CONTINUE TOP = 0 BOT = 0 L = K DO 50 K = 1,L TOP = TOP + I(K)*(I(K)+2)*E(K) 50 BOT = BOT+E(K) F = G2B2NK/12.D0/T * TOP/BOT RETURN END C OTHER functions C Trimer eq. S1=S2=S3=1/2, J12=J23=J31=J C Ref: Mossoba M. M. et al. J. Am. Chem. Soc.(1980) 102,6864. SUBROUTINE TRIMER(A,IPAR,T,F) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT DIMENSION A(1) IPAR = IPAR+1 AJ3=3.D0*A(IPAR)*HCBIK/T IF (AJ3.GT.EXPLMT) WRITE(IO,*) 1 'Exponent Warning in TRIMER: AJ*3 = ',AJ3 X=DEXP(AJ3) F=G2B2NK/(4.D0*T)*(1.D0+5.D0*X)/(1.D0+X) RETURN END C Hatfield's four-spin coupling in a cube (S1=S2=S3=S4=1/2) C Ref: Hall, J. W. et al. Inorg. Chem. (1977) 16, 1572. SUBROUTINE HATFLD(A,IPAR,T,F) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT DIMENSION A(1) IPAR = IPAR+1 AJ1=A(IPAR)*HCBIK/T IPAR = IPAR+1 AJ2 = A(IPAR)*HCBIK/T IPAR = IPAR+1 AJ3 = A(IPAR)*HCBIK/T IPAR = IPAR+1 AJ4 = A(IPAR)*HCBIK/T AJ1PJ2 = AJ1+AJ2 AJ1MJ2 = AJ1-AJ2 AJ3PJ4 = AJ3+AJ4 AJ3MJ4 = AJ3-AJ4 TT1=AJ1PJ2+2.D0*AJ3PJ4 Q= -(TT1)/2 T1=-(AJ1PJ2-2.D0*AJ3PJ4)/2 TT2=2.D0*DSQRT(AJ1MJ2*AJ1MJ2+AJ3MJ4*AJ3MJ4) T2=(AJ1PJ2+TT2)/2.D0 T3=(AJ1PJ2-TT2)/2.D0 TT3=2.D0*DSQRT((AJ1PJ2-AJ3PJ4)**2 + 3.D0*AJ3MJ4*AJ3MJ4) S1=(TT1+TT3)/2.D0 S2=(TT1-TT3)/2.D0 IF (DABS(Q).GT.EXPLMT) 1 WRITE(IO,*)'Exponent Warning in HATFLD--> Q =',Q Q = DEXP(-Q) T1 = DEXP(-T1) T2 = DEXP(-T2) T3 = DEXP(-T3) S1 = DEXP(-S1) S2 = DEXP(-S2) TOP = 5.D0*Q+T1+T2+T3 BOT = 5.D0*Q+3.D0*(T1+T2+T3)+S1+S2 CON = 2.D0*G2B2NK/T F = CON * TOP/BOT RETURN END C O. Kahn's eqn. for chain of dimers C Ref: Lorett, F. et al. Inorg. Chem. (1993) 32, 27. SUBROUTINE KAHN(A,IPAR,T,F,ITERM) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT DIMENSION A(1) IPAR = IPAR+1 AJ=A(IPAR) IPAR = IPAR+1 AJ2 = A(IPAR)*HCBIK/T CALL DIMER (S(ITERM),Z(ITERM),AJ,T,F) SSP1 = F/G2B2NK*3.D0*T X=AJ2*SSP1 U=1.D0/DTANH(X)-1.D0/X F=F*(1.D0+U)/(1.D0-U) RETURN END C Initial definitions for OTHER equations SUBROUTINE OTHPAR(NF,I) IMPLICIT REAL*8 (A-H,O-Z) COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK5/ IN, IO COMMON /BLK9/ NOTH 1 WRITE(IO,5) 5 FORMAT(' Select one of the following from the Other functions'/ 1 ' 1. User function'/ 1 ' 2. S=1/2 Three-spin system'/ 2 ' 3. Hatfield''s S=1/2 Four-spin system'/ 3 ' 4. Kahn''s expression') WRITE(IO,*) ' ' WRITE(IO,*) 'Enter Option number: ' READ(IN,*) NOTH GOTO (10,20,30,40), NOTH GOTO 1 10 WRITE(IO,*)'User defined function' CALL USRPAR(NF,I) GOTO 200 20 WRITE(IO,*) 'Three Spin S1=S2=S3= 1/2 system with single J:' CALL GETPAR('J-Trimr') GOTO 200 30 WRITE(IO,*) 'Hatfield''s 4 spin (S=1/2) system:' CALL GETPAR('J1-Tetra') CALL GETPAR('J2-Tetra') CALL GETPAR('J3-Tetra') CALL GETPAR('J4-Tetra') GOTO 200 40 WRITE(IO,*) 'Kahn''s expression for chain of dimers:' WRITE(IO,*)'Enter S1, S2 of dimer' READ(IN,*) S(I),Z(I) CALL GETPAR('J-Intra ') CALL GETPAR('J-Inter ') 200 RETURN END C Function selection for OTHER equations. C ITERM is necessary if NF(ITERM) value is used in any C of the functions used in this routine. SUBROUTINE OTHER(A,IPAR,T,F,ITERM) IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(1) COMMON /BLK5/ IN, IO COMMON /BLK6/ ITYPE, NFN, NF(10), NPAR COMMON /BLK9/ NOTH GOTO (10,20,30,40), NOTH 10 CALL USER(A,IPAR,T,F) GOTO 100 20 CALL TRIMER(A,IPAR,T,F) F=F/3 GOTO 100 30 CALL HATFLD(A,IPAR,T,F) F=F/4 GOTO 100 40 CALL KAHN(A,IPAR,T,F,ITERM) F=F/2 100 RETURN END C Both the subroutines USRPAR and USER need to be C modified to change user defined function. C USRPAR defines the names and initial values of minimization C parameters. This is done by CALL GETPAR statement once for C each parameter. C C Spin 1-1 dimer with isotropic D. Define parameters. SUBROUTINE USRPAR(NF,I) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK5/ IN, IO C C----------Modify the following portion of the routine to change ------ C----------to a new equation------------------------------------------- C WRITE(IO,*) 'S1=1, S2=1 dimer eq. including D' CALL GETPAR('J-1-1 ') CALL GETPAR('D-1-1 ') C C---------------------------------------------------------------------- RETURN END C Spin 1-1 dimer with isotropic D. Function evaluation SUBROUTINE USER(A,IPAR,T,F) IMPLICIT REAL*8 (A-H, O-Z) COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK5/ IN, IO COMMON /LMTS/ EXPLMT, NPRLMT, NFNLMT DIMENSION A(1) C---------------------------------------------------------------------- C This routine defines user specified function which may be C treated as a term in the susceptibility expression. C A - Parameter array. C IPAR - is the last parameter number used. If this routine C defines any new parameters, IPAR should be accordingly C increased before RETURNing to calling routine. C T - Temperature (K) C F - Susceptibility calculated by this routine C HCBIK = hc/k; B2NBIK = beta*beta*N/k C G2B2NK = g*g*beta*beta*N/k C HCNB2 = hc/(N*beta*beta) C HCNGB2 = hc/(N*g*beta*beta) C BK = beta*k BN = beta*N GBK=g*BK GBN=g*BN C If g is a minimization parameter, g = A(1), else g = G. C The condition G<1.E-6is used to check whether G is a minimzn. parm. C IN,IO = Input, Output unit numbers C EXPLMT = Maximum exponent value EXPLMT = DLOG(1.0D307) C---------------------------------------------------------------------- C E.g. In this routine, 2 parameters AJ and D are used to calculate C susceptibility (Xm) value. For each parameter used, IPAR is C incremented once. The final Xm value is defined as the variable F C which is returned. T is temperature supplied by the calling C routine. C---------------------------------------------------------------------- C-----The following portion is necessary to keep IPAR count consistent--- C---------------------------------------------------------------------- C IPAR = IPAR+1 AJ=A(IPAR) IPAR = IPAR+1 D = A(IPAR) C C---------------------------------------------------------------------- C The rest of the part defines the value of F (cald. suscep.) C---------------------------------------------------------------------- X = AJ*HCBIK/T D = D*HCBIK/T CON = G2B2NK/T/3 SJ = 6.D0*X TJ = 2.D0*X IF (DABS(SJ).GT.EXPLMT) 1 WRITE(IO,*)'Exponent Warning in USER--> SJ =',SJ E1 = DEXP(TJ - 2.D0*D/3.D0) E2 = DEXP(TJ + D/3.D0) E3 = DEXP(SJ - 2.D0*D) E4 = DEXP(SJ - D) E5 = DEXP(SJ+2.D0*D) TOP = 2.D0*E1+4.D0*E2+6.D0*E3+12.D0*(E4+E5) BOT = 1.D0+E1+E3+2.D0*(E2+E4+E5) C---------------------------------------------------------------------- C F is returned to calling routine C---------------------------------------------------------------------- F = CON * TOP/BOT/2 RETURN END C EQN writes the expression on file SUSCEP.OUT SUBROUTINE EQN IMPLICIT REAL*8 (A-H,O-Z) CHARACTER P*3,PP(10)*3 COMMON /BLK2/ HCBIK, B2NBIK, G, G2B2NK, HCNB2, HCNGB2 COMMON /BLK3/ S(10), Z(10), NS(10) COMMON /BLK6/ ITYPE, NFN, NF(10), NPAR COMMON /BLK9/ NOTH DATA PP/'A. ','B. ','C. ','D. ','E. ','F. ' 1 ,'G. ','H. ','I. ','J. '/ WRITE(2,*)'Equation type used for fitting of data:' IF(ITYPE.EQ.1) WRITE(2,*)' A + B + ...' IF(ITYPE.EQ.2) WRITE(2,*)' x(A) + (1-x)(B) + C + ...' IF(ITYPE.EQ.3) WRITE(2,*) 1 ' y [ x(A) + (1-x)(B) ] + (1-y) (C) + D + ...' IF(ITYPE.EQ.1)WRITE(2,*)'Where,' IF(ITYPE.EQ.2)WRITE(2,*)'Where, x is coefficient and' IF(ITYPE.EQ.3)WRITE(2,*)'Where, x and y are coefficients and' DO 100 I=1,NFN P=PP(I) N = NF(I) IF(N.EQ.1)WRITE(2,5)P,'Curie term',S(I) 5 FORMAT(1X,2A,' S=',F3.1) IF(N.EQ.2)WRITE(2,5)P,'Zero field term of',S(I) IF(N.EQ.3)WRITE(2,*)P,'Ising term' IF(N.EQ.4)WRITE(2,*)P,'Bleaney-Bower term' IF(N.EQ.5)WRITE(2,*)P,'Magnetization term' IF(N.EQ.6)WRITE(2,*)P,'Bonner Fisher term' IF(N.EQ.7)WRITE(2,*)P,'Modified Bonner Fisher term' IF(N.EQ.8)WRITE(2,20)P,'Intradimer term',S(I),Z(I) 20 FORMAT(1X,2A,' S1=',F3.1,' S2=',F3.1) 30 IF(N.NE.9)GOTO 40 IF(NOTH.EQ.1)WRITE(2,*)P,'User defined term' IF(NOTH.EQ.2)WRITE(2,*)P,'S=1/2 trimer term' IF(NOTH.EQ.3)WRITE(2,*)P,'Hatfield''s S=1/2 tetramer term' IF(NOTH.EQ.4)WRITE(2,20)P,'Kahn''s dimer chain term',S(I),Z(I) IF(NOTH.EQ.5)WRITE(2,*)P,'Heucher''s term' 40 IF(N.EQ.10) WRITE(2,*)P,'DIM term' IF(NS(I).NE.1) GOTO 100 IF(N.EQ.1) WRITE(2,*) 'Weiss correction (theta) applied' IF(N.GT.1) WRITE(2,*) 'Neighbor correction applied' 100 CONTINUE IF (DABS(G) .GT. 1.D-6) WRITE(2,*)'g-value fixed' IF (.NOT.(DABS(G) .GT. 1.D-6)) WRITE(2,*)'g-value varied' RETURN END C Set-up X-coordinate values for plot SUBROUTINE XDATA(T,NP,TC,NC,NPLOP,IPLOP) IMPLICIT REAL*8 (A-H, O-Z) REAL*4 XOBS,XCALC,YOBS,YCALC, XMIN,XMAX,YMIN,YMAX DIMENSION T(1),TC(1) COMMON /BLK7/ XOBS(300), XCALC(300), YOBS(300), YCALC(300) COMMON /BLK8/ XMIN, XMAX, YMIN, YMAX NC=300 XMIN = T(1) XMAX = T(1) DO 10 I=2,NP IF(XMAX .LT. T(I)) XMAX=T(I) IF(XMIN .GT. T(I)) XMIN=T(I) 10 CONTINUE IF(IPLOP.EQ.2) GOTO 40 XINT=(XMAX-XMIN)/(NC-1) DO 30 I=1,NC 30 TC(I)=(I-1)*XINT+XMIN GOTO 60 40 NC=NP DO 50 I=1,NC 50 TC(I)=T(I) 60 IF(NPLOP.EQ.5) GOTO 100 DO 70 I=1,NP 70 XOBS(I)=T(I) DO 80 I=1,NC 80 XCALC(I)=TC(I) GOTO 150 100 DO 110 I=1,NP 110 XOBS(I) = 1.D0/T(I) DO 120 I=1,NC 120 XCALC(I)=1.D0/TC(I) TEMP=XMIN XMIN=1.D0/XMAX XMAX=1.D0/TEMP 150 RETURN END C Set-up y-coordinate values for plot SUBROUTINE YDATA(T,XI,NP,TC,XIC,NC,NPLOP) IMPLICIT REAL*8 (A-H, O-Z) REAL*4 XOBS,XCALC,YOBS,YCALC, XMIN,XMAX,YMIN,YMAX COMMON /BLK7/ XOBS(300), XCALC(300), YOBS(300), YCALC(300) COMMON /BLK8/ XMIN, XMAX, YMIN, YMAX DIMENSION T(1), XI(1), TC(1), XIC(1) GOTO (10,30,50,70,10), NPLOP 10 DO 20 I=1,NP 20 YOBS(I)=XI(I) DO 25 I=1,NC 25 YCALC(I)=XIC(I) GOTO 100 30 T7=7.997D0 DO 40 I=1,NP 40 YOBS(I) = DSQRT(T7*XI(I)*T(I)) DO 45 I=1,NC 45 YCALC(I) = DSQRT(T7*XIC(I)*TC(I)) GOTO 100 50 DO 60 I=1,NP 60 YOBS(I) = XI(I)*T(I) DO 65 I=1,NC 65 YCALC(I)=XIC(I)*TC(I) GOTO 100 70 DO 80 I=1,NP 80 YOBS(I)=1.D0/XI(I) DO 85 I=1,NC 85 YCALC(I)=1.D0/XIC(I) 100 YMIN=YOBS(1) YMAX=YOBS(1) DO 110 I=2,NP IF(YMAX.LT.YOBS(I)) YMAX=YOBS(I) IF(YMIN.GT.YOBS(I)) YMIN=YOBS(I) 110 CONTINUE DO 120 I=1,NC IF(YMAX.LT.YCALC(I)) YMAX=YCALC(I) IF(YMIN.GT.YCALC(I)) YMIN=YCALC(I) 120 CONTINUE RETURN END C calling routine for Graphics plot SUBROUTINE SPLOT(NP,NC) COMMON /BLK7/ XOBS(300), XCALC(300), YOBS(300), YCALC(300) COMMON /BLK8/ XMIN, XMAX, YMIN, YMAX COMMON /BLK5/ IN, IO XDIV=(XMAX-XMIN)/5. YDIV=(YMAX-YMIN)/5. CALL XYPLOT(XOBS,YOBS,NP,XCALC,YCALC,NC,1, 1 XMIN,XMAX,XDIV,YMIN,YMAX,YDIV) RETURN END C Graphics plot routine for IBMPC compatibles using C Microsoft Fortran V 5.0 SUBROUTINE XYPLOT(XO,YO,M,XC,YC,N,L1, 1 XMIN,XMAX,XDIV,YMIN,YMAX,YDIV) C XC - X-COORDINATE ARRAY CALCULATED C YC - Y-COORDINATE ARRAY CALCULATED C XO - X-COORDINATE ARRAY EXPERIMENTAL C YO - Y-COORDINATE ARRAY EXPERIMENTAL C N - NUMBER OF CALCULATED POINTS C M - NUMBER OF EXPTAL DATA POINTS C L1,L2 - LINE TYPES FOR YO & YC CURVES C LINE TYPES : 1 - Solid 2 - Dotted 3 - Dashed C 4 - Dash-dot 5 - long dash 6- double dot C (X0,Y0) IS ORIGIN. (XL,YL) DIMENSIONS IN PLOTTER COORDS. C EXPECTS THE FONT FILE ROMAN.FON BE PRESENT IN THE CURRENT DIRECTORY. DIMENSION XC(1),YC(1),XO(1),YO(1) C ------------------------------------------------------- INCLUDE 'FGRAPH.FD' C -------------------------------------------------------- INTEGER*2 FX,FY,DUMMY,LS(6),X1,Y1,X2,Y2,MX,MY,NMX,NMY CHARACTER STR*5, STR1*2 RECORD /XYCOORD/ S RECORD /VIDEOCONFIG/ SCR FX(X)= (X-XMIN)*XS+X0 FY(Y)= NMY-((Y-YMIN)*YS+Y0) DATA LS/#FFFF, #AAAA, #CCCC, #EBEB, #7878, #8888/ C Initialize graphics mode CALL GETVIDEOCONFIG(SCR) NAD=SCR.ADAPTER C CGA HIGH RESOLUTION BW IF (NAD.EQ.#0002)MODE = 6 C HERCULES IF (NAD.EQ.#0020)MODE = 8 C EGA MONO & COLOR IF ((NAD.EQ.#0004 .OR. NAD.EQ.#0044) 1 .AND. SCR.MONITOR .EQ. #0001)MODE = 15 IF ((NAD.EQ.#0004 .OR. NAD.EQ.#0044) 1 .AND. SCR.MONITOR .NE. #0001)MODE = 16 C VGA, MCGA IF (NAD.EQ.#0008.OR.NAD.EQ.#0048.OR.NAD.EQ.#0010)MODE=17 DUMMY= SETVIDEOMODE(INT2(MODE)) C Define the plot dimensions CALL GETVIDEOCONFIG(SCR) NMX=SCR.NUMXPIXELS-1 NMY=SCR.NUMYPIXELS-1 X0=NMX*.16 Y0=NMY*.20 XL=NMX*.80 YL=NMY*.75 C X and Y Scale factors XS=XL/(XMAX-XMIN) YS=YL/(YMAX-YMIN) C Draw the first curve IF(L1.LT.1 .OR. L1.GT.6 )L1 =1 CALL SETLINESTYLE(LS(L1)) CALL MOVETO(FX(XC(1)),FY(YC(1)),S) DO 20 I=2,N 20 DUMMY=LINETO(FX(XC(I)),FY(YC(I))) C Draw the second curve (draw x marks) XCR = 2 DO 30 I=1,M CALL MOVETO(FX(XO(I))-XCR,FY(YO(I))-XCR,S) DUMMY=LINETO(FX(XO(I))+XCR,FY(YO(I))+XCR) CALL MOVETO(FX(XO(I))-XCR,FY(YO(I))+XCR,S) 30 DUMMY=LINETO(FX(XO(I))+XCR,FY(YO(I))-XCR) C Draw a box NH=12 NW=10 DUMMY=REGISTERFONTS('*.fon') DUMMY=SETFONT("t'ROMAN' h12 w10") MX=5*NW/2 MY=NH/2 40 CALL SETLINESTYLE(LS(1)) CALL MOVETO(FX(XMIN),FY(YMIN),S) DUMMY=LINETO(FX(XMIN),FY(YMAX)) DUMMY=LINETO(FX(XMAX),FY(YMAX)) DUMMY=LINETO(FX(XMAX),FY(YMIN)) DUMMY=LINETO(FX(XMIN),FY(YMIN)) C Draw ticks and values on X-axis L=ALOG10(XMAX-XMIN) ST=XMIN Y1=NMY-(Y0+NMY*.025) Y2=NMY-(Y0-NMY*.025) 50 CALL MOVETO(FX(ST),FY(YMIN),S) DUMMY=LINETO(FX(ST),Y1) A=ST/10.0**L WRITE(STR,'(F5.2)') A X2=FX(ST)-MX IF(X2.GT.INT2(NMX*.95)) X2=NMX*.95 CALL MOVETO(X2,Y2,S) CALL OUTGTEXT(STR) ST=ST+XDIV IF (ST.LE.XMAX) GOTO 50 IF (L.EQ.0) GOTO 55 CALL MOVETO(INT2(X0+XL-10*NW),INT2(Y0+YL+2.5*NH),S) CALL OUTGTEXT('X10') CALL MOVETO(INT2(X0+XL-7*NW),INT2(Y0+YL+2.0*NH),S) IF (L.GT.0 .AND. L.LE.9) WRITE(STR1,'(I1)')L IF (L.LT.0 .OR. L.GT.9) WRITE(STR1,'(I2)')L 54 CALL OUTGTEXT(STR1) C Draw ticks and values on Y-axis 55 L=ALOG10(YMAX-YMIN) ST=YMIN X1=X0+.008*NMX X2=X0-5.5*NW 60 CALL MOVETO(FX(XMIN),FY(ST),S) DUMMY=LINETO(X1,FY(ST)) A=ST/10.0**L WRITE(STR,'(F5.2)')A Y2=FY(ST)-MY IF (Y2.LT.INT2(0)) Y2=0 CALL MOVETO(X2,Y2,S) CALL OUTGTEXT(STR) ST=ST+YDIV IF (ST.LE.YMAX) GOTO 60 IF (L.EQ.0) GOTO 70 CALL MOVETO(INT2(X0-9*NW),INT2(NMY-(Y0+YL-1.5*NH)),S) CALL OUTGTEXT('X10') IF(L.GT.9 .OR. L.LT.0) WRITE(STR1,'(I2)')L IF(L.GT.0 .AND. L.LT.9) WRITE(STR1,'(I1)')L CALL MOVETO(INT2(X0-6*NW),INT2(NMY-(Y0+YL-NH)),S) CALL OUTGTEXT(STR1) C Wait for user's instruction to return 70 READ(*,*) CALL UNREGISTERFONTS() DUMMY=SETVIDEOMODE(INT2(-1)) RETURN END C Printer plot routine SUBROUTINE PPLOT(NP) CHARACTER*1 A(75),BLANK,STAR,CALC COMMON /BLK7/ X(300), XCALC(300), Y(300), Z(300) COMMON /BLK8/ XMIN, XMAX, YMIN, YMAX COMMON /BLK5/ IN, IO DATA BLANK, STAR, CALC /' ','*','c'/ NX = 74 NY =20 XS = NX/(XMAX-XMIN) YS = NY/(YMAX-YMIN) YINT = 1.D0/YS T2 = YMAX T1 = T2 NX1 = NX + 1 115 T1 = T1-YINT DO 120 J = 1,NX1 120 A(J) = BLANK NX1 = 0 DO 140 I=1,NP IF(Y(I) .LT.T1 .OR. Y(I) .GT. T2) GOTO 130 J = (X(I)-XMIN)* XS + 1.5 IF (J.GT.NX1) NX1=J A(J) = STAR 130 IF(Z(I).LT. T1 .OR. Z(I) .GT. T2) GOTO 140 J = (X(I)-XMIN)*XS + 1.5 IF (J.GT.NX1) NX1=J A(J) = CALC 140 CONTINUE WRITE(IO,150) (A(J),J=1,NX1) 150 FORMAT(4X,75A1) T2 = T1 IF (T1 .GT. YMIN) GOTO 115 READ(IN,160) A(1) 160 FORMAT(A) 200 RETURN END C Generates SUSCEP.HP file for plotting on HP7475A plotter. SUBROUTINE HPPLOT(NP,NC,NPLOP) INTEGER X0,Y0,XX,YY, XSHFT,YSHFT REAL INCH CHARACTER CH3*1 COMMON /BLK7/ XOBS(300), XCALC(300), YOBS(300), YCALC(300) COMMON /BLK8/ XMIN, XMAX, YMIN, YMAX COMMON /BLK5/ IN, IO IX(X)= (X-XMN)*XS+X0 IY(Y)= (Y-YMN)*YS+Y0 CH3 = CHAR(3) INCH=402 * 2.54 X0= 2.0 * INCH Y0= 1.5 * INCH XL= 8.5 * INCH YL= 6.0 * INCH XD=(XMAX-XMIN)/10. YD=(YMAX-YMIN)/10. WRITE(IO,*) 'The values from data are:' WRITE(IO,5) ' X-min =',XMIN,' X-max =',XMAX,' X-divn =',XD 5 FORMAT(3(A,F10.4)) WRITE(IO,*) 'Enter the values for plot' READ(IN,*) XMN, XMX, XD WRITE(IO,5) ' Y-min =',YMIN,' Y-max =',YMAX,' Y-divn =',YD WRITE(IO,*) 'Enter the values for plot' READ(IN,*) YMN, YMX, YD XS=XL/(XMX-XMN) YS=YL/(YMX-YMN) C Draw the curve OPEN(UNIT=3, FILE='SUSCEP.HP', STATUS='UNKNOWN') WRITE(IO,*)' HPPLOT file name is SUSCEP.HP' WRITE(3,10) C Initialize the plotter 10 FORMAT(1X,'IN;SP1;PU;') C Plot calculated curve DO 20 I=1,NC 20 WRITE(3,21)IX(XCALC(I)),IY(YCALC(I)) 21 FORMAT(1X,'PA',I5,',',I5,'PD') WRITE(3,22) 22 FORMAT(1X,'PU;') C Plot triangles for experimental points DO 25 I=1,NP 25 WRITE(3,26) IX(XOBS(I)),IY(YOBS(I)) 26 FORMAT(1X,'PA',I5,',',I5,'CI60,10WG60,0,360,10;') C Draw the box WRITE(3,21) IX(XMN),IY(YMN) WRITE(3,21) IX(XMX),IY(YMN) WRITE(3,21) IX(XMX),IY(YMX) WRITE(3,21) IX(XMN),IY(YMX) WRITE(3,21) IX(XMN),IY(YMN) WRITE(3,22) C Write X-ticks WRITE(3,27) 27 FORMAT(1X,'SI0.2,0.3;') XTB=INT(XMN/XD)*XD IF(XTB .LT. XMN) XTB=XTB+XD NTX=(XMX-XTB)/XD +1 DO 50 I=1,NTX XX=(XTB-XMN)*XS+X0 WRITE(3,30) XX,Y0 30 FORMAT(1X,'PA',I5,',',I5,'XT;') IXTB=XTB+0.5 IF(NPLOP.EQ.5) GOTO 40 IF(IXTB.LT.10) WRITE(3,31) IXTB,CH3 IF(IXTB.GE.10 .AND. IXTB.LT.100) WRITE(3,32) IXTB,CH3 IF(IXTB.GE.100) WRITE(3,33) IXTB,CH3 31 FORMAT(1X,'CP-.5,-1.5;LB',I1,A) 32 FORMAT(1X,'CP-1.,-1.5;LB',I2,A) 33 FORMAT(1X,'CP-1.5,-1.5;LB',I3,A) GOTO 50 40 WRITE(3,41) XTB,CH3 41 FORMAT(1X,'CP-1.5,-1.5;LB',F3.2,A) 50 XTB = XTB + XD C Write Y-ticks YTB=INT(YMN/YD)*YD IF(YTB.LT.YMN) YTB = YTB+YD NTY=(YMX-YTB)/YD+1 DO 100 I=1,NTY YY=(YTB-YMN)*YS+Y0 WRITE(3,60) X0,YY 60 FORMAT(1X,'PA',I5,',',I5,'YT;') GOTO (80,70,70,90,80), NPLOP 70 WRITE(3,71) YTB,CH3 71 FORMAT(1X,'CP-3.5,-.5,LB',F3.1,A) GOTO 100 80 IYTB=YTB*1000 WRITE(3,81) IYTB,CH3 81 FORMAT(1X,'CP-3.5,-.5,LB',I3,A) GOTO 100 90 IYEXP=0 IF (YMX.LE.0) GOTO 91 IYEXP=LOG10(YMX) 91 YEXP = 10.**IYEXP YTBE=YTB/YEXP WRITE(3,71) YTBE,CH3 100 YTB=YTB+YD C X-Axis title XSHFT=X0+XL/2. WRITE(3,110) XSHFT,Y0 110 FORMAT(1X,'PA',I5,',',I5,';') IF(NPLOP.LT.5) WRITE(3,111) CH3 IF(NPLOP.EQ.5) WRITE(3,112) CH3, CH3, CH3 111 FORMAT(1X,'CP-3,-3;LBT (K) ',A) 112 FORMAT(1X,'CP-4,-3;LB1/T (K',A,'CP0,.5;SI.18,.25;LB-1',A, 1 'SI.2,.3; CP0,-.5;LB)',A) WRITE(3,195) C Y-Axis title YSHFT=Y0+YL/2 WRITE(3,130) X0, YSHFT 130 FORMAT(1X,'PA',I5,',',I5,'SI.2,.3;DI0,1;') GOTO (140,150,160,170,140), NPLOP 140 WRITE(3,141) 141 FORMAT(1X,'CP-6,-3;') WRITE(3,180) WRITE(3,142) CH3,CH3,CH3,CH3,CH3,CH3,CH3 142 FORMAT(1X,'LB (x10',A,'CP0,.5;SI.18,.25;LB3',A, 1 'SI.2,.3;CP0,-.5;LBcm',A,'CP0,.5;SI.18,.25;LB3',A,/ 2 ' SI.2,.3;CP0,-.5;LBmol',A,'CP0,.5;SI.18,.25;LB-1', 3 A,'SI.2,.3;CP0,-.5LB) ',A) GOTO 200 150 WRITE(3,151) 151 FORMAT(1X,'CP-4,-3;') C MU WRITE(3,185) WRITE(3,152) CH3 152 FORMAT(1X,'LB (B.M.) ',A) GOTO 200 160 WRITE(3,161) 161 FORMAT(1X,'CP-10,-3;') WRITE(3,180) WRITE(3,162) CH3,CH3,CH3,CH3,CH3,CH3,CH3 162 FORMAT(1X,'LB T (x10',A,'CP0,.5;SI.18,.25;LB3',A, 1 'SI.2,.3;CP0,-.5;LBcm',A,'CP0,.5;'/' SI.18,.25;LB3',A, 2 ' SI.2,.3;CP0,-.5;'/' LBmol',A,'CP0,.5;SI.18,.25;LB-1', 3 A,'SI.2,.3;CP0,-.5LBK) ',A) GOTO 200 170 WRITE(3,171) 171 FORMAT(1X,'CP-8,-3;') WRITE(3,180) WRITE(3,172) CH3,IYEXP,CH3,CH3,CH3,CH3 172 FORMAT(1X,'CP-3,0LB1/ (x10',A,'CP0,.5;SI.18,.25;LB',I1,A, 1 'SI.2,.3;CP0,-.5;LB cm',A,'CP0,.5;'/' SI.18,.25;LB-3',A, 2 ' SI.2,.3;CP0,-.5;LBmol)',A) GOTO 200 C Write Chi 180 FORMAT(1X,'PR-1206,0;'/ 1 ' PR-213,24PD;PR-1,2PR1,4PR4,4PR6,-2PR3,-4PR1,-6;',/ 1 ' PR-1,-6PR-4,-4PR-7,0PR-9,4PR-4,8PR-1,12PR1,12PR4,12;',/ 2 ' PR9,16PR15,16PR78,44PR4,8PR0,4PR-1,4PR-3,-8PR-3,0;',/ 3 ' PR-6,4PR-1,4PR-1,8PR1,8PR7,4PR6,-2PR6,-6PR1,-8;',/ 4 ' PR1,-12PR0,-8PR-3,-12PR-87,-52PR-12,-12PR-9,-16;',/ 5 ' PR-3,-8PR0,-8PR3,-4PR4,0PU;', 6 ' PR103,4PD;PR-9,-12PR-108,144PR9,12PR108,-144PU;') C Write mu 185 FORMAT(1X,'PR-1206,0;'/ 1 ' PR0,0PD;PR-120,40,0,20,42,-8,12,4,4,4,1,12;',/ 1 ' PR-1,8,-4,8,-7,8,-46,24,0,20,51,-20,4,2,1,6;'/ 2 ' PR-1,8,-4,8,4,0,4,-2,4,-6,1,-8,-1,-8,-4,-8;',/ 3 ' PR-6,-6,6,-6,7,-12,1,-16,0,-8,-3,-12,-4,-4;',/ 4 ' PR-4,-4,-24,6,90,-40,0,-20PU;') C Up arrow 190 FORMAT(1X,'PR-1206,0;','PR-51,0PD;PR0,120PR18,-20PR-27,60'/ 1 'PR-27,-60PR18,20PR0,-120PR18,0PU;') C Right arrow 195 FORMAT(1X,'PR100,50PD;PR90,0PR-15,-24PR45,36;'/ 1 ' PR-45,36PR15,-24PR-90,0PR0,-24PU;') 200 WRITE(3,210) 210 FORMAT(1X,'PU;SP0;') CLOSE(UNIT=3,STATUS='KEEP') RETURN END C C Ref: Error Estimation using the Sequential Simplex Method in C Nonlinear Least Squares Data Analysis C Gregory R. Phillips and Edward M. Eyring C Anal. Chem. 1988, vol 60, pp 738-741 C and references therein SUBROUTINE ERREST(XMIN, SIGMA,N,NP,IER) C C XMIN is converged parameter set by simplex method C P final simplex; N - No of parameters C NP - No. of observed points C IER = 1 if the routine fails IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION XMIN(12),MIN(12),SIGMA(12) DIMENSION B(12,12),BI(12,12),Q(12,12),A(12) 1 ,PBAR(12),Y(13),Y0I(12),BIA(12),BBA(12) 2 ,PP(20,21),YY(13) COMMON /ERRCAL/ P(20,21) DATA ZERO,HALF,TWO,THR,FOUR/0.0D0,0.5D0,2.0D0,3.0D0,4.0D0/ C DO 10 I=1,N 10 SIGMA(I)=ZERO IF(NP.LE.N) GOTO 550 IER=0 ROUND=1.D-12 NN=N+1 NM1=N-1 DNN=FLOAT(NN) C DO 30 I=1,NN DO 20 J=1,N 20 MIN(J)=P(J,I) 30 Y(I)=FN(MIN) C 100 ROUND=ROUND*10.D0 IF(ROUND.GT.1.0D-5) GOTO 550 ITER=0 IEXPN=1 C CALCULATE V0 110 ILO=1 YLO=Y(1) DO 200 I=1,NN IF(Y(I).GE.YLO) GOTO 200 YLO=Y(I) ILO=I 200 CONTINUE IF(IEXPN.EQ.0) GOTO 310 C CHECK WHETHER EXPANSION IS NECESSARY DO 220 I=1,NN IF(I.EQ.ILO) GOTO 220 IF((Y(I)-YLO)/YLO .LT. ROUND) GOTO 230 220 CONTINUE GOTO 280 C EXPANSION 230 DO 250 I=1,N Z=ZERO DO 240 J=1,NN 240 Z=Z+P(I,J) 250 PBAR(I)=Z/NN DO 270 I=1,NN DO 260 J=1,N P(J,I)=PBAR(J)+TWO*(P(J,I)-PBAR(J)) 260 MIN(J)=P(J,I) 270 Y(I)=FN(MIN) GOTO 110 C SAVE THE SIMPLEX 280 DO 290 I=1,N YY(I)=Y(I) DO 290 J=1,N 290 PP(J,I)=P(J,I) C CALCULATION BEGINS 300 IEXPN=0 GOTO 110 C MOVE COLUMN ILO TO COLUMN NN 310 DO 320 J=1,N 320 MIN(J)=P(J,ILO) Y0=Y(ILO) DO 330 I=ILO,N I1=I+1 Y(I)=Y(I1) DO 330 J=1,N 330 P(J,I)=P(J,I1) DO 340 J=1,N 340 P(J,NN)=MIN(J) Y(NN)=Y0 C CONSTRUCT A, B, Q MATRICES DO 360 I=1,N DO 350 J=1,N Q(J,I)=(P(J,I)-P(J,NN)) 350 MIN(J)=(P(J,I)+P(J,NN))*HALF Y0I(I)=FN(MIN) A(I)=FOUR*Y0I(I)-Y(I)-THR*Y0 360 B(I,I)=TWO*(Y(I)+Y0-TWO*Y0I(I)) DO 380 I=1,NM1 I1=I+1 DO 380 J=I1,N DO 370 K=1,N 370 MIN(K)=(P(K,I)+P(K,J))*HALF YIJ=FN(MIN) B(I,J)=TWO*(YIJ+Y0-Y0I(I)-Y0I(J)) 380 B(J,I)=B(I,J) CALL INVERT(B,BI,N,EPS) IF (EPS.LT.ZERO) GOTO 100 DO 400 I=1,N BIA(I)=ZERO DO 400 K=1,N 400 BIA(I)=BIA(I)+BI(I,K)*A(K) DO 410 I=1,N BBA(I)=ZERO DO 410 K=1,N 410 BBA(I)=BBA(I)+BI(I,K)*BIA(K) ABBA=ZERO DO 420 K=1,N 420 ABBA=ABBA+A(K)*BBA(K) IF(ABBA.GT. 0.25D0) ITER=ITER+1 IF(ITER.LE.4) GOTO 450 C TOO MANY ITERATIONS. EXPAND SIMPLEX. DO 430 I=1,NN Y(I)=YY(I) DO 430 J=1,N 430 P(J,I)=PP(J,I) GOTO 100 C CORRECT THETA0 450 DO 470 J=1,N QBAJ=ZERO DO 460 K=1,N 460 QBAJ=QBAJ+Q(J,K)*BIA(K) P(J,NN)=P(J,NN)-QBAJ*HALF 470 MIN(J)=P(J,NN) Y(NN)=FN(MIN) IF (ABBA.GT. 0.25D0) GOTO 300 C CALCULATE SIGMA 500 S2=Y(NN)/(NP-N) DO 520 I=1,N DO 520 J=1,N B(I,J)=ZERO DO 520 K=1,N 520 B(I,J)=B(I,J)+BI(I,K)*Q(J,K) DO 540 I=1,N SIGMA(I)=ZERO DO 530 K=1,N 530 SIGMA(I)=SIGMA(I)+Q(I,K)*B(K,I) IF(SIGMA(I).LT.ZERO) GOTO 100 540 SIGMA(I)=DSQRT(SIGMA(I)*S2) DO 545 I=1,N 545 XMIN(I)=MIN(I) RETURN 550 IER=1 RETURN END C Calls the SSP subroutine GJRD for inversion. C Convergence criterian = (Largest element of B) * 1.e-6 SUBROUTINE INVERT(B,BI,N,EPS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION B(12,12),BI(12,12),TEMP(144) IO=6 EPS=0.D0 K=0 DO 5 J=1,N DO 5 I=1,N IF (EPS.LT.DABS(B(I,J))) EPS=DABS(B(I,J)) K=K+1 5 TEMP(K)=B(I,J) EPS=EPS*1.D-6 CALL GJRD(TEMP,N,EPS) C IF (EPS.EQ.-1.0D0) WRITE(IO,*)' Error in matrix inversion' K=0 DO 10 J=1,N DO 10 I=1,N K=K+1 10 BI(I,J)=TEMP(K) RETURN END C IBM 1130 - SCIENTIFIC SUBROUTINE PACKAGE SUBROUTINE GJRD (A,N,EPS) C The following lines are modified for double precision C all print statements are suppressed C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(1) ,L(50), M(50) ABS(X)=DABS(X) C--------------------------------------------------------------- C N IS THE ORDER OF THE MATRIX A(I,K) TO BE INVERTED ,EPS IS A C TOLERANCE FOR ACCEPTANCE OF THE SINGULARITY OF THE GIVEN MATRIX C IF MATRIX SINGULAR THEN EPS=-1 C DETERMINATION OF THE PIVOT ELEMENT C IPR=6 IF(N-50) 1,1,99 1 CONTINUE C SEARCH FOR LARGEST ELEMENT D=1.0 NK=-N DO 80 K=1,N NK=NK+N L(K)=K M(K)=K KK=NK+K BIGA=A(KK) DO 20 J=K,N IZ=N*(J-1) DO 20 I=K,N IJ=IZ+I 10 IF( ABS(BIGA)- ABS(A(IJ))) 15,20,20 15 BIGA=A(IJ) L(K)=I M(K)=J 20 CONTINUE C INTERCHANGE ROWS J=L(K) IF(J-K) 35,35,25 25 KI=K-N DO 30 I=1,N KI=KI+N HOLD=-A(KI) JI=KI-K+J A(KI)=A(JI) 30 A(JI) =HOLD C INTERCHANGE COLUMNS 35 I=M(K) IF(I-K) 45,45,38 38 JP=N*(I-1) DO 40 J=1,N JK=NK+J JI=JP+J HOLD=-A(JK) A(JK)=A(JI) 40 A(JI) =HOLD C DIVIDE COLUMN BY MINUS PIVOT (VALUE OF PIVOT ELEMENT IS C CONTAINED IN BIGA) 45 IF(ABS(BIGA)- EPS) 46,46,48 46 EPS = - 1.0 C WRITE (IPR,47) K,BIGA 47 FORMAT(24H0MATRIX SINGULAR* STAGE,I6,4X5HPIVOTE20.10) RETURN 48 DO 55 I=1,N IF(I-K) 50,55,50 50 IK=NK+I A(IK)=A(IK)/(-BIGA) 55 CONTINUE C REDUCE MATRIX DO 65 I=1,N IK=NK+I HOLD=A(IK) IJ=I-N DO 65 J=1,N IJ=IJ+N IF(I-K) 60,65,60 60 IF(J-K) 62,65,62 62 KJ=IJ-I+K A(IJ)=HOLD*A(KJ)+A(IJ) 65 CONTINUE C DIVIDE ROW BY PIVOT KJ=K-N DO 75 J=1,N KJ=KJ+N IF(J-K) 70,75,70 70 A(KJ)=A(KJ)/BIGA 75 CONTINUE C PRODUCT OF PIVOTS D=D*BIGA C REPLACE PIVOT BY RECIPROCAL A(KK)=1.0/BIGA 80 CONTINUE C FINAL ROW AND COLUMN INTERCHANGE K=N 100 K=(K-1) IF(K) 150,150,105 105 I=L(K) IF(I-K) 120,120,108 108 JQ=N*(K-1) JR=N*(I-1) DO 110 J=1,N JK=JQ+J HOLD=A(JK) JI=JR+J A(JK)=-A(JI) 110 A(JI) =HOLD 120 J=M(K) IF(J-K) 100,100,125 125 KI=K-N DO 130 I=1,N KI=KI+N HOLD=A(KI) JI=KI-K+J A(KI)=-A(JI) 130 A(JI) =HOLD GO TO 100 150 RETURN 99 CONTINUE C WRITE(IPR,101)N 101 FORMAT(39H0***** SUBROUTINE GJRA ' N = 50 N =,I10,/) STOP END C C C ALGORITHM AS 47 APPLIED STATISTICS(J.R.STATIST.SOC C), C (1971) VOL 20,NO.3 C SUBROUTINE NELMIN(N,START,MIN,YNEWLO,REQMIN,STEP,KONGVE,ICOUNT) DOUBLE PRECISION START(N),MIN(N),YNEWLO,REQMIN,STEP(N), 1 P(20,21),PSTAR(20),P2STAR(20),PBAR(20),Y(20), 2 DN,DNN,Z,SUM,SUMM,YLO,RCOEFF,YSTAR,ECOEFF,Y2STAR,CCOEFF, 3 CURMIN,DEL,FN COMMON /ERRCAL/ P C DATA RCOEFF/1.D0/,ECOEFF/2.D0/,CCOEFF/5.D-1/ C REFLECTION,EXTENSION AND CONTRACTION COEFFICIENTS. C C VALIDITY CHECKS ON INPUT PARAMETERS. KCOUNT=ICOUNT ICOUNT=0 IF(REQMIN.LE.0.D0)ICOUNT=ICOUNT-1 IF(N.GT.20)ICOUNT=ICOUNT-10 IF(KONGVE.LE.0)ICOUNT=ICOUNT-100 IF(ICOUNT.LT.0)RETURN C JCOUNT=KONGVE DN=FLOAT(N) NN=N+1 DNN=FLOAT(NN) DEL=1.D0 C C CONSTRUCTION OF INITIAL SIMPLEX C 1001 DO 1 I=1,N 1 P(I,NN)=START(I) Z=FN(START) Y(NN)=Z SUM=Z SUMM=Z*Z DO 2 J=1,N START(J)=START(J)+STEP(J)*DEL DO 3 I=1,N 3 P(I,J)=START(I) Z=FN(START) Y(J)=Z SUM=SUM+Z SUMM=SUMM+Z*Z 2 START(J)=START(J)-STEP(J)*DEL C C SIMPLEX CONSTRUCTION COMPLETE C C FIND HIGHEST AND LOWEST Y VALUES.YNEWLO(=Y(IHI)) INDICATES C THE VERTEX OF THE SIMPLEX TO BE REPLACED C 1000 YLO=Y(1) YNEWLO=YLO ILO=1 IHI=1 DO 5 I=2,NN IF(Y(I).GE.YLO)GO TO 4 YLO=Y(I) ILO=I 4 IF(Y(I).LE.YNEWLO)GO TO 5 YNEWLO=Y(I) IHI=I 5 CONTINUE SUM=SUM-YNEWLO SUMM=SUMM-YNEWLO*YNEWLO C C CALCULATE PBAR,THE CENTROID OF THE SIMPLEX VERTICES C EXCEPTING THAT WITH Y VALUE YNEWLO DO 7 I=1,N Z=0.D0 DO 6 J=1,NN 6 Z=Z+P(I,J) Z=Z-P(I,IHI) 7 PBAR(I)=Z/DN C C REFLECTION THRO' THE CENTROID C DO 8 I=1,N 8 PSTAR(I)=(1.D0+RCOEFF)*PBAR(I)-RCOEFF*P(I,IHI) YSTAR=FN(PSTAR) ICOUNT=ICOUNT+1 IF(YSTAR.GE.YLO) GO TO 12 C C SUCCESSFUL REFLECTION,SO EXTENSION C DO 9 I=1,N 9 P2STAR(I)=ECOEFF*PSTAR(I)+(1.D0-ECOEFF)*PBAR(I) Y2STAR=FN(P2STAR) ICOUNT=ICOUNT+1 C C RETAIN EXTENSION OR CONTRACTION C IF(Y2STAR.GE.YLO)GO TO 19 10 DO 11 I=1,N 11 P(I,IHI)=P2STAR(I) Y(IHI)=Y2STAR GO TO 900 C C NO EXTENSION C 12 L=0 DO 13 I=1,NN IF(Y(I).GT.YSTAR)L=L+1 13 CONTINUE IF(L.GT.1)GO TO 19 IF(L.EQ.0)GO TO 15 C C CONTRACTION OR REFLECTION SIDE OF THE CENTROID C DO 14 I=1,N 14 P(I,IHI)=PSTAR(I) Y(IHI)=YSTAR C C CONTRACTION ON THE Y(IHI) SIDE OF THE CENTROID C 15 DO 16 I=1,N 16 P2STAR(I)=CCOEFF*P(I,IHI)+(1.D0-CCOEFF)*PBAR(I) Y2STAR=FN(P2STAR) ICOUNT=ICOUNT+1 IF(Y2STAR.LE.Y(IHI))GO TO 10 C C CONTRACT WHOLE SIMPLEX C SUM=0.D0 SUMM=0.D0 DO 18 J=1,NN DO 17 I=1,N P(I,J)=(P(I,J)+P(I,ILO))*0.5D0 17 MIN(I)=P(I,J) Y(J)=FN(MIN) SUM=SUM+Y(J) 18 SUMM=SUMM+Y(J)*Y(J) ICOUNT=ICOUNT+NN GO TO 901 C C RETAIN REFLECTION C 19 DO 20 I=1,N 20 P(I,IHI)=PSTAR(I) Y(IHI)=YSTAR 900 SUM=SUM+Y(IHI) SUMM=SUMM+Y(IHI)*Y(IHI) 901 JCOUNT=JCOUNT-1 IF(JCOUNT.NE.0) GO TO 1000 C C CHECK TO SEE IF MINIMUM IS REACHED C IF(ICOUNT.GT.KCOUNT)GO TO 22 JCOUNT=KONGVE CURMIN=(SUMM-(SUM*SUM)/DNN)/DN C C CURMIN IS THE VARIANCE OF N+1 FN VALUES AT THE VERTICES C IF(CURMIN.GE.REQMIN)GO TO 1000 C C FACTORIAL TESTS TO CHECK THAT YNEWLO IS A LOCAL MINIMUM C 22 DO 23 I=1,N 23 MIN(I)=P(I,IHI) YNEWLO=Y(IHI) IF(ICOUNT.GT.KCOUNT)GOTO 30 DO 24 I=1,N DEL=STEP(I)*1.D-3 MIN(I)=MIN(I)+DEL Z=FN(MIN) IF(Z.LT.YNEWLO)GO TO 25 MIN(I)=MIN(I)-DEL-DEL Z=FN(MIN) IF(Z.LT.YNEWLO)GO TO 25 24 MIN(I)=MIN(I)+DEL C ADDED PART TO ESTIMATE THE VARIATION IN VARIABLES 30 DO 32 I=1,N START(I)=0.D0 DO 32 J=1,NN IF (J.EQ.IHI) GOTO 32 Z=DABS(P(I,J)-P(I,IHI)) IF(START(I).LT.Z) START(I)=Z 32 CONTINUE C THE MAXIMUM DIFFERENCES IN PARAMETERS ARE IN START(I) RETURN C C RESTART PROCEDURE C 25 DO 26 I=1,N 26 START(I)=MIN(I) DEL=1.D-3 ICOUNT=ICOUNT+10000 GO TO 1001 END