C@(#)=================================================================== C@(#) C@(#) C@(#) FILE = tranfit.f C@(#) C@(#) --------------- VERSION = 1.15 C@(#) | SCCS FILE | C@(#) | SUMMARY | CURRENT CHECKOUT DATE = 08/10/94 C@(#) --------------- at 18:17:09 C@(#) DATE OF NEWEST DELTA = 08/10/94 C@(#) at 18:17:09 C@(#) SCCS file name = /users/chemkin/SCCS/s.tranfit.f C@(#)=================================================================== C23456789012345678901234567890123456789012345678901234567890123456789012 C/////////////////////////////////////////////////////////////////////// C SUBROUTINE TRANFT C C WRITTEN BY: C ROBERT J. KEE C COMPUTATIONAL MECHANICS DIVISION C SANDIA NATIONAL LABORATORIES C LIVERMORE, CA 94550 C (415) 294-3272 C C////////////////////////////////////////////////////////////////////// C C VERSION 3.3 C C CHANGES FROM PREVIOUS VERSION: C 1. Restructured data in OMEG12 C 2. Changed REAL*8 to DOUBLE PRECISION C 3. Changed POLFIT and PCOEF to call single and double precision C SLATEC subroutines POLFIT,DPOLFT and PCOEF,DPCOEF C 4. Change vms open statements C CHANGES FROM VERSION 1.3 C 1. Change THIGH to 3500 C 2. Change name to TRANFT to conform to ANSI standard C 3. Add "unix" and "snla" change blocks C CHANGES FROM VERSION 1.4 C 1. modify OPEN statements for unix C CHANGES FOR VERSION 1.6 C 1. Find THIGH and TLOW from species thermodynamic data C CHANGES FOR VERSION 1.9 C 1. Change binary file to include version, precision, error C status, and required length of work arrays C 2. Allow user input from a file. C 3. Implement non-formatted transport data input C CHANGES FOR VERSION 3.0 (3/15/94 F. Rupley) C 1. DOS/PC compatibility effort includes adding file names to C OPEN statements, removing unused variables in CALL lists, C unusued but possibly initialized variables. C CHANGES FOR VERSION 3.1 (6/9/94 H. Moffat) C 1. Changed Avrog Number to 1986 CODATA recommendations. C 2. PI and variables based on PI C is now accurate up to DP machine prec. C CHANGES FOR VERSION 3.2 (8/10/94 H.K. Moffat) C 1. Got rid of LEXIST variable that was uninitialized. C 2. 'tran.inp' is no longer openned - there is no need for it. C CHANGES FOR VERSION 3.3 (1/19/95 F. Rupley) C 1. Add integer error flag to CKINIT call list. C 2. Appended "polfit.f", as tranfit is the only code which C uses it. C///////////////////////////////////////////////////////////////////// C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C PARAMETER (LINKCK=10, LINKTP=35, LTRAN=31, LOUT=6, 1 KDIM=99, NFDIM=165, NO=4, NT=50, NCHECK=1000, 2 LENICK=10000, LENRCK=10000, LENCCK=100, MAXFIT=7, 3 NRANGE=2, MAXTMP=3) C DIMENSION EPS(KDIM), SIG(KDIM), DIP(KDIM), POL(KDIM), ZROT(KDIM), 1 NLIN(KDIM), WT(KDIM), CV(KDIM), FITWT(NT), FITRES(NT), 2 ALOGT(NT), XLA(NT), XETA(NT), XD(NT), FITWRK(NFDIM), 3 COFLAM(NO,KDIM), COFETA(NO,KDIM), COFD(NO,KDIM,KDIM), 4 COFTD(NO,KDIM,10), ICKWRK(LENICK), RCKWRK(LENRCK), 5 VALUE(6), KTDIF(KDIM), NTEMP(KDIM), TEMP(MAXTMP,KDIM), 6 THERM(MAXFIT, NRANGE, KDIM) C CHARACTER CCKWRK(LENCCK)*16, KSYM(KDIM)*16, LINE*80, VERS*16, 1 PREC*16 LOGICAL IERR, KERR COMMON /TPPARS/ PI, RU, PATM, BOLTZ, EPSIL, FDTCGS, FATCM, 1 DIPMIN C DATA NLIN/KDIM*NCHECK/, EMAXL,EMAXE,EMAXD,EMAXTD/4*0.0/, 1 KERR/.FALSE./ C C #ifdef DOUBLE_PRECISION PI = DACOS(-1.0D0) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION PI = ACOS(-1.0) #endif /* SINGLE_PRECISION */ EPSIL = 0.0 FDTCGS = 1.0D-18 FATCM = 1.0D+08 DIPMIN = 1.0D-20 C OPEN (LINKCK, FORM='UNFORMATTED', STATUS='OLD', 1 FILE='chem.bin') OPEN (LTRAN, FORM='FORMATTED', STATUS='OLD', 1 FILE='tran.dat') OPEN (LINKTP, FORM='UNFORMATTED', STATUS='UNKNOWN', 1 FILE='tran.bin') C C OPEN (LDATA,FORM='UNFORMATTED',STATUS='UNKNOWN', C 1 FILE='tran.inp') C C WRITE VERSION NUMBER C VERS = '3.3' WRITE (LOUT, 15) VERS(:3) 15 FORMAT( 1/' TRANFT: Transport property fitting,', 2/' CHEMKIN-II Version ',A,', January 1995', #ifdef DOUBLE_PRECISION 3/' DOUBLE PRECISION') PREC = 'DOUBLE' #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION 3/' SINGLE PRECISION') PREC = 'SINGLE' #endif /* SINGLE_PRECISION */ C WRITE (LOUT, 8310) C C INITIALIZE CHEMKIN C CALL CKINIT (LENICK, LENRCK, LENCCK, LINKCK, LOUT, ICKWRK, 1 RCKWRK, CCKWRK, IFLAG) CLOSE (LINKCK) IF (IFLAG .GT. 0) STOP CALL CKINDX (ICKWRK, RCKWRK, MM, KK, II, NFIT) IF (KK .GT. KDIM) THEN WRITE (LOUT, 8300) KDIM STOP ENDIF C CALL CKSYMS (CCKWRK, LOUT, KSYM, IERR) KERR = KERR.OR.IERR CALL CKWT (ICKWRK, RCKWRK, WT) CALL CKRP (ICKWRK, RCKWRK, RU, RUC, PATM) P = PATM C C Value of Avrogadro's number is from 1986 CODATA C recommended values. C (6.0221367(36)E23 mol-1) C #ifdef DOUBLE_PRECISION BOLTZ = RU / 6.0221367D23 #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION BOLTZ = RU / 6.0221367E23 #endif /* SINGLE_PRECISION */ C CALL CKATHM (MAXFIT, NRANGE, ICKWRK, RCKWRK, MAXTMP, 1 NTEMP, TEMP, THERM) TLOW = TEMP(1,1) THIGH = TEMP(NTEMP(1),1) DO 10 K = 2, KK THIGH = MIN (THIGH, TEMP (NTEMP(K), K)) TLOW = MAX (TLOW, TEMP (1,K)) 10 CONTINUE DT = (THIGH-TLOW) / (NT-1) C WRITE (LOUT, 20) TLOW, THIGH 20 FORMAT (/,' DATA HAS BEEN FIT OVER THE RANGE: TLOW=',F9.2, 1 ', THIGH=',F9.2) C C READ THE TRANSPORT DATA BASE C WRITE (LOUT, 7020) LIN = LTRAN C 50 CONTINUE READ (LIN, '(A)', END=500) LINE ILEN = IPPLEN(LINE) IF (ILEN .GT. 0) THEN K1 = IFIRCH(LINE(:ILEN)) K2 = K1 + INDEX(LINE(K1:ILEN),' ') - 1 CALL CKCOMP (LINE(K1:K2-1), KSYM, KK, KFIND) C IF (KFIND .GT. 0) THEN CALL CKXNUM (LINE(K2:ILEN), 6, LOUT, NVAL, VALUE, IERR) KERR = KERR.OR.IERR IF (IERR) WRITE (LOUT, 8000) WRITE (LOUT, 7010) LINE C NLIN(KFIND) = INT(VALUE(1)) EPS(KFIND) = VALUE(2) SIG(KFIND) = VALUE(3) DIP(KFIND) = VALUE(4) POL(KFIND) = VALUE(5) ZROT(KFIND) = VALUE(6) ENDIF ENDIF GO TO 50 C 500 CONTINUE C C IF (LIN .EQ. LTRAN) THEN C LIN = LDATA C GO TO 50 C ENDIF C DO 600 K = 1, KK IF (NLIN(K) .EQ. NCHECK) THEN DO 750 J = K, KK IF (NLIN(J) .EQ. NCHECK) WRITE (LOUT, 8010) KSYM(J) 750 CONTINUE KERR = .TRUE. ENDIF 600 CONTINUE C IF (.NOT. KERR) THEN C C FIT THE CONDUCTIVITIES AND VISCOSITIES C CALL LAMFIT (KK, NT, NO, LOUT, WT, SIG, EPS, DIP, ZROT, NLIN, 1 P, TLOW, DT, ALOGT, FITRES, FITWT, FITWRK, XLA, 2 XETA, CV, ICKWRK, RCKWRK, COFLAM, COFETA, EMAXL, 3 EMAXE) C C FIT THE DIFFUSION COEFFICIENTS C CALL DIFFIT (KK, NT, NO, KDIM, LOUT, WT, SIG, EPS, DIP, POL, 1 P, TLOW, DT, ALOGT, FITRES, FITWT, FITWRK, XD, 2 COFD, EMAXD) C C FIT THE THERMAL DIFFUSION RATIOS C CALL THMFIT (KK, NT, NO, KDIM, LOUT, WT, SIG, EPS, DIP, POL, 1 TLOW, DT, ALOGT, FITRES, FITWT, FITWRK, XD, 2 NLITE, KTDIF, COFTD, EMAXTD) C C PRINT THE FITS C WRITE (LOUT, 7035) 1 ' COEFFICIENTS FOR SPECIES CONDUCTIVITIES',EMAXL WRITE (LOUT, 8200) WRITE (LOUT, 8100) (KSYM(K), (COFLAM(N,K),N=1,NO), K=1,KK) C WRITE (LOUT, 7035) 1 ' COEFFICIENTS FOR SPECIES VISCOSITIES',EMAXE WRITE (LOUT, 8200) WRITE (LOUT, 8100) (KSYM(K), (COFETA(N,K),N=1,NO), K=1,KK) C WRITE (LOUT, 7035) 1 ' COEFFICIENTS FOR SPECIES DIFFUSION COEFFICIENTS',EMAXD DO 2300 J = 1, KK WRITE (LOUT, 8200) WRITE (LOUT, 8110) 1 (KSYM(J), KSYM(K), (COFD(N,J,K),N=1,NO), K=1,J) 2300 CONTINUE C WRITE (LOUT, 7035) 1 ' COEFFICIENTS FOR THERMAL DIFFUSION RATIOS',EMAXTD DO 2400 M = 1, NLITE K = KTDIF(M) WRITE (LOUT, 8200) WRITE (LOUT, 8110) 1 (KSYM(K), KSYM(J), (COFTD(N,J,M),N=1,NO), J=1,KK) 2400 CONTINUE ELSE WRITE (LOUT, '(/A)') 1 ' WARNING...THERE IS AN ERROR IN THE TRANSPORT LINKING FILE' ENDIF C C WRITE THE LINKING FILE C WRITE (LINKTP) VERS, PREC, KERR C LENIMC = 4*KK + NLITE LENRMC = (19 + 2*NO + NO*NLITE)*KK + (15+NO)*KK**2 C WRITE (LINKTP) LENIMC, LENRMC, NO, KK, NLITE WRITE (LINKTP) PATM, (WT(K), EPS(K), SIG(K), DIP(K), 1 POL(K), ZROT(K), NLIN(K), K=1,KK), 2 ((COFLAM(N,K),N=1,NO),K=1,KK), 4 ((COFETA(N,K),N=1,NO),K=1,KK), 5 (((COFD(N,J,K),N=1,NO),J=1,KK),K=1,KK), 6 (KTDIF(N),N=1,NLITE), 7 (((COFTD(N,J,L),N=1,NO),J=1,KK),L=1,NLITE) CLOSE (LINKTP) C RETURN C C FORMATS C 7000 FORMAT (A) 7010 FORMAT (1X,A) 7020 FORMAT (///,' TRANSPORT PARAMETERS FROM DATA BASE:',/) 7035 FORMAT (///A/' MAXIMUM FITTING ERROR = ', 1PE12.3) 8000 FORMAT (10X, 'ERROR IN CKXNUM READING FROM TRANSPORT DATA BASE') 8010 FORMAT (10X, 'ERROR...TRANSPORT DATA NOT GIVEN FOR ',A10) 8100 FORMAT (2X, A10, 4E12.3) 8110 FORMAT (2X, A10, A10, 4E12.3) 8200 FORMAT (1X, ' ') 8300 FORMAT (10X, 'ERROR...THE TRANSPORT FITTING CODE IS DIMENSIONED 1FOR ONLY', I3, ' SPECIES') 8310 FORMAT(///,' OUTPUT FROM THE TRANSPORT PROPERTY FITTING PACKAGE' 1,//) END C C----------------------------------------------------------------------- C SUBROUTINE LAMFIT (KK, NT, NO, NOUT, WT, SIG, EPS, DIP, ZROT,NLIN, 1 P, TLOW, DT, ALOGT, FITRES, FITWT, FITWRK, XLA, 2 XETA, CV, ICKWRK, RCKWRK, COFLAM, COFETA,EMAXL, 3 EMAXE) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C DIMENSION WT(*), SIG(*), EPS(*), DIP(*), ZROT(*), NLIN(*), CV(*), 1 FITRES(*), FITWRK(*), ALOGT(*), FITWT(*), XLA(*), 2 XETA(*), ICKWRK(*), RCKWRK(*), COFLAM(NO,*), 3 COFETA(NO,*) C COMMON /TPPARS/ PI, RU, PATM, BOLTZ, EPSIL, FDTCGS, FATCM, DIPMIN C FCON = 0.5 * FDTCGS**2 * FATCM**3 / BOLTZ C PFAC = PATM / P DO 1000 K = 1, KK DST = FCON * DIP(K)**2 / (EPS(K) * SIG(K)**3) HELPE = 2.6693D-5 * SQRT(WT(K)) / SIG(K)**2 HELPD = 2.6280D-3 * PFAC / (SQRT(WT(K)) * SIG(K)**2) C DO 300 N = 1, NT T = TLOW + (N-1)*DT PRUT = P / (RU * T) TR = T / EPS(K) CALL CKCVML (T, ICKWRK, RCKWRK, CV) ALOGT(N) = LOG(T) DII = SQRT(T**3)*HELPD / OMEG12(1,TR,DST) XETA(N) = SQRT(T)*HELPE / OMEG12(2,TR,DST) RODET = DII * WT(K) * PRUT / XETA(N) AA = 2.5 - RODET CALL PARKER (T, EPS(K), PARK) IF (NLIN(K) .EQ. 2) THEN BB = ZROT(K)*PARK + 2.0*(5.0/2.0 + RODET)/PI CROCT = 1.0 ELSE BB = ZROT(K)*PARK + 2.0*(5.0/3.0 + RODET)/PI CROCT = 2.0/3.0 ENDIF FTRA = 2.5 * (1.0 - 2.0 * CROCT * AA / (BB*PI)) FROT = RODET * (1.0 + 2.0 * AA / (BB*PI)) FVIB = RODET IF (NLIN(K) .EQ. 0) THEN FAKTOR = 5.0/2.0 * 1.5*RU ELSEIF (NLIN(K) .EQ. 1) THEN FAKTOR = (FTRA*1.5 + FROT)*RU + FVIB*(CV(K)-2.5*RU) ELSEIF (NLIN(K) .EQ. 2) THEN FAKTOR = (FTRA+FROT)*1.5*RU + FVIB*(CV(K)-3.0*RU) ENDIF XLA(N) = LOG( XETA(N)/WT(K) * FAKTOR) XETA(N) = LOG(XETA(N)) 300 CONTINUE C C FIT CONDUCTIVITY C FITWT(1) = -1.0 EPSL = EPSIL C #ifdef DOUBLE_PRECISION CALL DPOLFT (NT, ALOGT, XLA, FITWT, NO-1, NORD, EPSL, FITRES, 1 IERR, FITWRK) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL POLFIT (NT, ALOGT, XLA, FITWT, NO-1, NORD, EPSL, FITRES, 1 IERR, FITWRK) #endif /* SINGLE_PRECISION */ C EMAXL = MAX (EMAXL, EPSL) IF (IERR .NE. 1) THEN WRITE (NOUT,7000) STOP ENDIF C CCC = 0.0 #ifdef DOUBLE_PRECISION CALL DPCOEF (NORD, CCC, COFLAM(1,K), FITWRK) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL PCOEF (NORD, CCC, COFLAM(1,K), FITWRK) #endif /* SINGLE_PRECISION */ C C FIT VISCOSITY C FITWT(1) = -1.0 EPSL = EPSIL #ifdef DOUBLE_PRECISION CALL DPOLFT #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL POLFIT #endif /* SINGLE_PRECISION */ C 1 (NT, ALOGT, XETA, FITWT, NO-1, NORD, EPSL, FITRES, IERR, FITWRK) C EMAXE = MAX (EMAXE, EPSL) IF (IERR .NE. 1) THEN WRITE (NOUT,7000) STOP ENDIF C CCC = 0.0 #ifdef DOUBLE_PRECISION CALL DPCOEF (NORD, CCC, COFETA(1,K), FITWRK) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL PCOEF (NORD, CCC, COFETA(1,K), FITWRK) #endif /* SINGLE_PRECISION */ C 1000 CONTINUE C 7000 FORMAT(9X,'ERROR IN POLFIT WHILE FITTING VISCOSITY OR ', 1'CONDUCTIVITY') RETURN END C SUBROUTINE DIFFIT (KK, NT, NO, KDIM, NOUT, WT, SIG, EPS, 1 DIP, POL, P, TLOW, DT, ALOGT, FITRES, FITWT, 2 FITWRK, XD, COFD, EMAXD) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C DIMENSION WT(*), SIG(*), EPS(*), DIP(*), POL(*), FITRES(*), 1 FITWT(*), FITWRK(*), ALOGT(*), XD(*), COFD(NO,KDIM,*) C COMMON /TPPARS/ PI, RU, PATM, BOLTZ, EPSIL, FDTCGS, FATCM, DIPMIN C DO 50 N = 1, NT T = TLOW + (N-1)*DT ALOGT(N) = LOG(T) 50 CONTINUE C FCON = 0.5 * FDTCGS**2 * FATCM**3 / BOLTZ C DO 1000 J = 1, KK DO 1000 K = 1, J SIGM = 0.5 * (SIG(J)+SIG(K)) EPSM = SQRT(EPS(J)*EPS(K)) IF (DIP(J).LT.DIPMIN .AND. DIP(K).GT.DIPMIN) THEN C C K IS POLAR, J IS NONPOLAR C DST = 0. XI = 1.0 + 1 POL(J) * FCON * DIP(K)**2 * SQRT(EPS(K)/EPS(J)) / 2 ( 2.0 * SIG(J)**3 * EPS(K) * SIG(K)**3) C ELSEIF (DIP(J).GT.DIPMIN .AND. DIP(K).LT.DIPMIN) THEN C C J IS POLAR, K IS NONPOLAR C DST = 0. XI = 1.0 + 1 POL(K) * FCON * DIP(J)**2 * SQRT(EPS(J)/EPS(K)) / 2 (2.0 * SIG(K)**3 * EPS(J) * SIG(J)**3) C C ELSE C C NORMAL CASE, EITHER BOTH POLAR OR BOTH NONPOLAR C DST = FCON * DIP(J) * DIP(K) / (EPSM * SIGM**3) XI = 1.0 ENDIF C SIGM = SIGM * XI**(-1.0/6.0) EPSM = EPSM * XI**2 HELP1 = (WT(J)+WT(K)) / (2.0*WT(J)*WT(K)) DO 500 N = 1, NT T = TLOW + (N-1)*DT TR = T / EPSM HELP2 = 0.002628 * SQRT(HELP1*T**3) XD(N) = HELP2*PATM / (OMEG12(1,TR,DST) * SIGM**2 * P) XD(N) = XD(N) * 1 D12(WT(J),WT(K),T,EPS(J),EPS(K),SIG(J),SIG(K),DST) XD(N) = LOG(XD(N)) 500 CONTINUE C FITWT(1) = -1.0 EPSL = EPSIL C #ifdef DOUBLE_PRECISION CALL DPOLFT #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL POLFIT #endif /* SINGLE_PRECISION */ 1 (NT, ALOGT, XD, FITWT, NO-1, NORD, EPSL, FITRES, IERR, FITWRK) C EMAXD = MAX (EMAXD, EPSL) IF (IERR .NE. 1) THEN WRITE (NOUT,7000) J,K STOP ENDIF C CCC = 0.0 C #ifdef DOUBLE_PRECISION CALL DPCOEF (NORD, CCC, COFD(1,K,J), FITWRK) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL PCOEF (NORD, CCC, COFD(1,K,J), FITWRK) #endif /* SINGLE_PRECISION */ 1000 CONTINUE C C FILL OUT THE LOWER TRIANGLE C DO 2000 K = 1, KK-1 DO 2000 J = K+1, KK DO 2000 N = 1, NO COFD(N,J,K) = COFD(N,K,J) 2000 CONTINUE RETURN 7000 FORMAT (10X, 'ERROR IN FITTING', I3, ',', I3, 1 'DIFFUSION COEFFICIENT') END C SUBROUTINE THMFIT (KK, NT, NO, KDIM, NOUT, WT, SIG, EPS, DIP, 1 POL, TLOW, DT, AT, FITRES, FITWT, FITWRK, TD, 2 NLITE, KTDIF, COFTD, EMAXTD) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C DIMENSION WT(*), SIG(*), EPS(*), DIP(*), POL(*), KTDIF(*), 1 COFTD(NO,KDIM,*), FITRES(*), FITWT(*), FITWRK(*), 2 AT(*), TD(*), FITAST(7), FITBST(7), FITCST(7) C COMMON /TPPARS/ PI, RU, PATM, BOLTZ, EPSIL, FDTCGS, FATCM, DIPMIN C DATA FITAST/ .1106910525D+01, -.7065517161D-02,-.1671975393D-01, 1 .1188708609D-01, .7569367323D-03,-.1313998345D-02, 2 .1720853282D-03/ C DATA FITBST/ .1199673577D+01, -.1140928763D+00,-.2147636665D-02, 1 .2512965407D-01, -.3030372973D-02,-.1445009039D-02, 2 .2492954809D-03/ C DATA FITCST/ .8386993788D+00, .4748325276D-01, .3250097527D-01, 1 -.1625859588D-01, -.2260153363D-02, .1844922811D-02, 2 -.2115417788D-03/ C NLITE = 0 DO 50 N = 1, NT T = TLOW + (N-1)*DT AT(N) = T 50 CONTINUE C DO 1100 J = 1, KK C IF (WT(J) .LE. 5.0) THEN NLITE = NLITE + 1 KTDIF(NLITE) = J EPSJ = EPS(J) * BOLTZ SIGJ = SIG(J) * 1.0D-8 POLJ = POL(J) * 1.0D-24 POLJST = POLJ / SIGJ**3 C DO 1000 K = 1, KK EPSK = EPS(K) * BOLTZ SIGK = SIG(K) * 1.0D-8 DIPK = DIP(K) * 1.0D-18 DIPKST = DIPK / SQRT(EPSK*SIGK**3) EKOEJ = EPSK / EPSJ TSE = 1.0 + 0.25*POLJST*DIPKST**2*SQRT(EKOEJ) EOK = TSE**2 * SQRT(EPS(J)*EPS(K)) WTKJ= (WT(K)-WT(J)) / (WT(K)+WT(J)) C DO 500 N = 1, NT TSLOG = LOG(AT(N) / EOK) CALL HORNER (7, FITAST, TSLOG, ASTAR) CALL HORNER (7, FITBST, TSLOG, BSTAR) CALL HORNER (7, FITCST, TSLOG, CSTAR) TD(N) = 7.5 * WTKJ * (2.0*ASTAR + 5.0) * 1 (6.0*CSTAR - 5.0)/ 2 (ASTAR * (16.0*ASTAR - 12.0*BSTAR + 55.0)) 500 CONTINUE C FITWT(1) = -1.0 EPSL = EPSIL C #ifdef DOUBLE_PRECISION CALL DPOLFT #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL POLFIT #endif /* SINGLE_PRECISION */ 1 (NT, AT, TD, FITWT, NO-1, NORD, EPSL, FITRES, IERR, FITWRK) C EMAXTD = MAX (EMAXTD, EPSL) IF (IERR .NE. 1) THEN WRITE (NOUT, 7000) J,K STOP ENDIF C CCC = 0.0 C #ifdef DOUBLE_PRECISION CALL DPCOEF (NORD, CCC, COFTD(1,K,NLITE), FITWRK) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION CALL PCOEF (NORD, CCC, COFTD(1,K,NLITE), FITWRK) #endif /* SINGLE_PRECISION */ C 1000 CONTINUE ENDIF 1100 CONTINUE 7000 FORMAT (10X, 'ERROR IN FITTING THE ', I3, ',', I3, 1 ' THERMAL DIFFUSION RATIO') RETURN END C SUBROUTINE HORNER (N, A, X, P) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C DIMENSION A(*) NM1 = N-1 B = A(N) DO 10 I = 1, NM1 B = A(N-I) + B*X 10 CONTINUE P = B RETURN END C #ifdef DOUBLE_PRECISION DOUBLE PRECISION FUNCTION D12 1 (W1, W2, T, EPS1, EPS2, SIG1, SIG2, DST) IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION REAL FUNCTION D12 (W1, W2, T, EPS1, EPS2, SIG1, SIG2, DST) IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C C CORRECTION OF D(1,2). C REFERENCE: MARRERO AND MASON,J.PHYS CHEM REF DAT.1,3(1972) C SUMW = W1+W2 SIG12 = 0.5 * (SIG1 + SIG2) TR11 = T / EPS1 TR22 = T / EPS2 TR12 = T / SQRT(EPS1*EPS2) C CALL OMEGXX (2, TR11, OM2F11, DST) CALL OMEGXX (2, TR22, OM2F22, DST) CALL OMEGXX (1, TR12, OM1F12, DST) CALL OMEGXX (2, TR12, OM2F12, DST) CALL OMEGXX (3, TR12, BST12, DST) CALL OMEGXX (4, TR12, CST12, DST) AST12 = OM2F12 / OM1F12 C H = (SIG1/SIG2)**2 * SQRT(2.0*W2/SUMW) * 2.0*W1**2/(SUMW*W2) P1 = H * OM2F11/OM1F12 C H = (SIG2/SIG1)**2 * SQRT(2.0*W1/SUMW) * 2.0*W2**2/(SUMW*W1) P2 = H * OM2F22/OM1F12 C P12 = 15.0 * ((W1-W2)/SUMW)**2 + 8.0*W1*W2*AST12/SUMW**2 C H = 2.0/(W2*SUMW) * SQRT(2.0*W2/SUMW) * OM2F11/OM1F12 * 1 (SIG1/SIG12)**2 Q1 = ((2.5-1.2*BST12)*W1**2 + 1 3.0*W2**2 + 1.6*W1*W2*AST12)*H C H = 2.0/(W1*SUMW) * SQRT(2.0*W1/SUMW) * OM2F22/OM1F12 * 1 (SIG2/SIG12)**2 Q2 = ((2.5-1.2*BST12)*W2**2 + 1 3.0*W1**2 + 1.6*W1*W2*AST12)*H C H = ((W1-W2)/SUMW)**2 * (2.5-1.2*BST12)*15.0 + 1 4.0*W1*W2*AST12/SUMW**2 * (11.0 - 2.4*BST12) Q12 = H + 1.6*SUMW*OM2F11*OM2F22 / (SQRT(W1*W2)*OM1F12**2) * 1 SIG1**2*SIG2**2/SIG12**4 D12 = ((6.0*CST12-5.0)**2/10.0) * (P1+P2+P12) / 1 (Q1+Q2+Q12) + 1.0 C RETURN END C SUBROUTINE PARKER (T, EPS, P) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C C TEMPERATURE DEPENDENCE OF THE ROTATIONAL COLLISION NUMBERS C C REF: BRAU, C.A., JONKMAN, R.M., "CLASSICAL THEORY C OF ROTATIONAL RELAXATION IN DIATOMIC GASES", C JCP,VOL52,P477(1970) C #ifdef DOUBLE_PRECISION DATA PI32O2 /2.7841639984158539D+00/, $ P2O4P2 /4.4674011002723397D+00/, $ PI32 /5.5683279968317078D+00/ #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION DATA PI32O2 /2.7841639984158539E+00/, $ P2O4P2 /4.4674011002723397E+00/, $ PI32 /5.5683279968317078E+00/ #endif /* SINGLE_PRECISION */ D = EPS / T DR = EPS / 298.0 P = (1.0 + PI32O2*SQRT(DR) + P2O4P2*DR + PI32*DR**1.5) / 1 (1.0 + PI32O2*SQRT(D) + P2O4P2*D + PI32*D**1.5) RETURN END C SUBROUTINE INTERP (ARG, X, Y, VAL) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C C//// QUADRATIC INTERPOLATION ////////////////////////////////////////// C C DIMENSION X(3),Y(3) VAL1 = Y(1) + (ARG-X(1))*(Y(2)-Y(1)) / (X(2)-X(1)) VAL2 = Y(2) + (ARG-X(2))*(Y(3)-Y(2)) / (X(3)-X(2)) FAC1 = (ARG-X(1)) / (X(2)-X(1)) / 2.0 FAC2 = (X(3)-ARG) / (X(3)-X(2)) / 2.0 IF (ARG .GE. X(2)) THEN VAL = (VAL1*FAC2+VAL2) / (1.0+FAC2) ELSE VAL = (VAL1+VAL2*FAC1) / (1.0+FAC1) ENDIF RETURN END C #ifdef DOUBLE_PRECISION DOUBLE PRECISION FUNCTION OMEG12 (N, TR, DR) IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION REAL FUNCTION OMEG12 (N, TR, DR) IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C C//// CALC. OF COLLISION INTEGRALS FOR A KRIEGER-12-6-3-POTENTIAL ////// C DIMENSION VERT(3), ARG(3), VAL(3), TSTERN(37), DELTA(8), O(37,8), 1 P(37,8) DATA TSTERN/.1,.2,.3,.4,.5,.6,.7,.8,.9,1.,1.2,1.4,1.6,1.8,2.,2.5, 1 3.,3.5,4.,5.,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.,25., 2 30.,35.,40.,50.,75.,100./ DATA DELTA/0.,.25,.5,.75,1.,1.5,2.,2.5/ DATA ((O(I,L),L=1,8),I=1,10) / 1 4.008 , 4.002 , 4.655 , 5.52 , 6.454 , 8.214 , 9.824 ,11.31 , 2 3.130 , 3.164 , 3.355 , 3.721 , 4.198 , 5.23 , 6.225 , 7.160 , 3 2.649 , 2.657 , 2.77 , 3.002 , 3.319 , 4.054 , 4.785 , 5.483 , 4 2.314 , 2.32 , 2.402 , 2.572 , 2.812 , 3.386 , 3.972 , 4.539 , 5 2.066 , 2.073 , 2.14 , 2.278 , 2.472 , 2.946 , 3.437 , 3.918 , 6 1.877 , 1.885 , 1.944 , 2.06 , 2.225 , 2.628 , 3.054 , 3.747 , 7 1.729 , 1.738 , 1.79 , 1.893 , 2.036 , 2.388 , 2.763 , 3.137 , 8 1.6122, 1.622 , 1.67 , 1.76 , 1.886 , 2.198 , 2.535 , 2.872 , 9 1.517 , 1.527 , 1.572 , 1.653 , 1.765 , 2.044 , 2.35 , 2.657 , * 1.44 , 1.45 , 1.49 , 1.564 , 1.665 , 1.917 , 2.196 , 2.4780/ DATA ((O(I,L),L=1,8),I=11,20) / 1 1.3204, 1.33 , 1.364 , 1.425 , 1.51 , 1.72 , 1.956 , 2.199 , 2 1.234 , 1.24 , 1.272 , 1.324 , 1.394 , 1.573 , 1.777 , 1.99 , 3 1.168 , 1.176 , 1.202 , 1.246 , 1.306 , 1.46 , 1.64 , 1.827 , 4 1.1166, 1.124 , 1.146 , 1.185 , 1.237 , 1.372 , 1.53 , 1.7 , 5 1.075 , 1.082 , 1.102 , 1.135 , 1.181 , 1.3 , 1.441 , 1.592 , 6 1.0006, 1.005 , 1.02 , 1.046 , 1.08 , 1.17 , 1.278 , 1.397 , 7 .95 , .9538, .9656, .9852, 1.012 , 1.082 , 1.168 , 1.265 , 8 .9131, .9162, .9256, .9413, .9626, 1.019 , 1.09 , 1.17 , 9 .8845, .8871, .8948, .9076, .9252, .972 , 1.03 , 1.098 , * .8428, .8446, .850 , .859 , .8716, .9053, .9483, .9984/ DATA ((O(I,L),L=1,8),I=21,30) / 1 .813 , .8142, .8183, .825 , .8344, .8598, .8927, .9316, 2 .7898, .791 , .794 , .7993, .8066, .8265, .8526, .8836, 3 .7711, .772 , .7745, .7788, .7846, .8007, .822 , .8474, 4 .7555, .7562, .7584, .7619, .7667, .78 , .7976, .8189, 5 .7422, .743 , .7446, .7475, .7515, .7627, .7776, .796 , 6 .72022, .7206, .722 , .7241, .7271, .7354, .7464, .76 , 7 .7025, .703 , .704 , .7055, .7078, .7142, .7228, .7334, 8 .68776, .688, .6888, .6901, .6919, .697 , .704 , .7125, 9 .6751, .6753, .676 , .677 , .6785, .6827, .6884, .6955, * .664 , .6642, .6648, .6657, .6669, .6704, .6752, .681 / DATA ((O(I,L),L=1,8),I=31,37) / 1 .6414, .6415, .6418, .6425, .6433, .6457, .649 , .653 , 2 .6235, .6236, .6239, .6243, .6249, .6267, .629 , .632 , 3 .60882, .6089, .6091, .6094, .61 , .6112, .613 , .6154, 4 .5964, .5964, .5966, .597 , .5972, .5983, .600 , .6017, 5 .5763, .5763, .5764, .5766, .5768, .5775, .5785, .58 , 6 .5415, .5415, .5416, .5416, .5418, .542 , .5424, .543 , 7 .518 , .518 , .5182, .5184, .5184, .5185, .5186, .5187/ C DATA ((P(I,L),L=1,8),I=1,10) / 1 4.1 , 4.266 , 4.833 , 5.742 , 6.729 , 8.624 ,10.34 ,11.890 , 2 3.263 , 3.305 , 3.516 , 3.914 , 4.433 , 5.57 , 6.637 , 7.618 , 3 2.84 , 2.836 , 2.936 , 3.168 , 3.511 , 4.329 , 5.126 , 5.874 , 4 2.531 , 2.522 , 2.586 , 2.749 , 3.004 , 3.64 , 4.282 , 4.895 , 5 2.284 , 2.277 , 2.329 , 2.46 , 2.665 , 3.187 , 3.727 , 4.249 , 6 2.084 , 2.081 , 2.13 , 2.243 , 2.417 , 2.862 , 3.329 , 3.786 , 7 1.922 , 1.924 , 1.97 , 2.072 , 2.225 , 2.641 , 3.028 , 3.435 , 8 1.7902, 1.795 , 1.84 , 1.934 , 2.07 , 2.417 , 2.788 , 3.156 , 9 1.682 , 1.689 , 1.733 , 1.82 , 1.944 , 2.258 , 2.596 , 2.933 , * 1.593 , 1.60 , 1.644 , 1.725 , 1.84 , 2.124 , 2.435 , 2.746 / DATA ((P(I,L),L=1,8),I=11,20) / 1 1.455 , 1.465 , 1.504 , 1.574 , 1.67 , 1.913 , 2.181 , 2.45 , 2 1.355 , 1.365 , 1.4 , 1.461 , 1.544 , 1.754 , 1.989 , 2.228 , 3 1.28 , 1.289 , 1.321 , 1.374 , 1.447 , 1.63 , 1.838 , 2.053 , 4 1.222 , 1.231 , 1.26 , 1.306 , 1.37 , 1.532 , 1.718 , 1.912 , 5 1.176 , 1.184 , 1.209 , 1.25 , 1.307 , 1.45 , 1.618 , 1.795 , 6 1.0933, 1.1 , 1.119 , 1.15 , 1.193 , 1.304 , 1.435 , 1.578 , 7 1.039 , 1.044 , 1.06 , 1.083 , 1.117 , 1.204 , 1.31 , 1.428 , 8 .9996, 1.004 , 1.016 , 1.035 , 1.062 , 1.133 , 1.22 , 1.32 , 9 .9699, .9732, .983 , .9991, 1.021 , 1.08 , 1.153 , 1.236 , * .9268, .9291, .936 , .9473, .9628, 1.005 , 1.058 , 1.12 / DATA ((P(I,L),L=1,8),I=21,30) / 1 .8962, .8979, .903 , .9114, .923 , .9545, .9955, 1.044 , 2 .8727, .8741, .878 , .8845, .8935, .918 , .9505, .9893, 3 .8538, .8549, .858 , .8632, .8703, .890 , .9164, .9482, 4 .8379, .8388, .8414, .8456, .8515, .868 , .8895, .916 , 5 .8243, .8251, .8273, .8308, .8356, .8493, .8676, .89 , 6 .8018, .8024, .8039, .8065, .810 , .820 , .8337, .8504, 7 .7836, .784 , .7852, .7872, .7899, .7976, .808 , .8212, 8 .7683, .7687, .7696, .771 , .7733, .7794, .788 , .7983, 9 .7552, .7554, .7562, .7575, .7592, .764 , .771 , .7797, * .7436, .7438, .7445, .7455, .747 , .7512, .757 , .7642/ DATA ((P(I,L),L=1,8),I=31,37) / 1 .71982, .72 , .7204, .7211, .7221, .725 , .7289, .7339, 2 .701 , .7011, .7014, .702 , .7026, .7047, .7076, .7112, 3 .68545, .6855, .686 , .686 , .6867, .6883, .6905, .693 , 4 .6723, .6724, .6726, .673 , .6733, .6745, .676 , .6784, 5 .651 , .651 , .6512, .6513, .6516, .6524, .6534, .6546, 6 .614 , .614 , .6143, .6145, .6147, .6148, .6148, .6147, 7 .5887, .5889, .5894, .59 , .5903, .5901, .5895, .5885/ C IF (DR.LT.-.00001 .OR. DR.GT.2.5 .OR. TR.LT..09 .OR. TR.GT.500. 1 .OR. (ABS(DR).GT.1.0D-5 .AND. TR.GT.75.0)) THEN WRITE (6,'(A)') 'COLLISION INTEGRAL UNDEFINED' OMEG12 = 0.0 RETURN ENDIF C IF (TR .GT. 75.0) THEN IF (N .EQ. 1) THEN OMEG12 = 0.623 - 0.136D-2*TR + 0.346D-5*TR*TR 1 -0.343D-8*TR*TR*TR ELSEIF (N .EQ. 2) THEN OMEG12 = 0.703 - 0.146D-2*TR + 0.357D-5*TR*TR 1 -0.343D-8*TR*TR*TR ENDIF RETURN ENDIF C IF (TR .LE. 0.2) THEN II = 2 ELSE II = 37 DO 30 I = 2, 36 IF (TSTERN(I-1).LT.TR .AND. TSTERN(I).GE.TR) II=I 30 CONTINUE ENDIF C IF (ABS(DR) .GE. 1.0D-5) THEN IF (DR .LE. 0.25) THEN KK = 2 ELSE KK = 7 DO 10 K = 2, 7 IF (DELTA(K-1).LT.DR .AND. DELTA(K) .GE. DR) KK=K 10 CONTINUE ENDIF DO 50 I = 1, 3 DO 40 K = 1, 3 ARG(K) = DELTA(KK-2+K) IF (N .EQ. 1) THEN VAL(K) = O(II-2+I, KK-2+K) ELSEIF (N .EQ. 2) THEN VAL(K) = P(II-2+I, KK-2+K) ENDIF 40 CONTINUE CALL INTERP (DR, ARG, VAL, VERT(I)) 50 CONTINUE DO 60 I = 1, 3 ARG(I) = TSTERN(II-2+I) 60 CONTINUE CALL INTERP (TR, ARG, VERT, OMEG12) C ELSE DO 80 I = 1, 3 ARG(I) = TSTERN(II-2+I) IF (N .EQ. 1) THEN VAL(I) = O(II-2+I, 1) ELSEIF (N .EQ. 2) THEN VAL(I) = P(II-2+I, 1) ENDIF 80 CONTINUE CALL INTERP (TR, ARG, VAL, OMEG12) ENDIF RETURN END C SUBROUTINE OMEGXX (N, TR, OM, DST) C #ifdef DOUBLE_PRECISION IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER (I-N) #endif /* DOUBLE_PRECISION */ #ifdef SINGLE_PRECISION IMPLICIT REAL (A-H, O-Z), INTEGER (I-N) #endif /* SINGLE_PRECISION */ C C//// CALCULATION OF THE REDUCED COLLISION INTEGRALS /////////////////// C IF (N.EQ.1 .OR. N.EQ.2) THEN OM = OMEG12 (N, TR, DST) ELSEIF (N .EQ. 3) THEN IF (TR .LE. 5.0) THEN OM = 1.36 - 0.223*TR + 0.0613*TR**2 -0.00554*TR**3 ELSE OM = 1.095 ENDIF ELSEIF (N .EQ. 4) THEN IF (TR .LE. 5.0) THEN OM =0.823 + 0.0128*TR + 0.0112*TR**2 -0.00193*TR**3 ELSE OM = .9483 ENDIF ELSE WRITE (6,'(10X,A)') 'OMEGXX IS CALLED WITH WRONG PARAMETERS' ENDIF RETURN END #ifdef DOUBLE_PRECISION SUBROUTINE DP1VLU (L, NDER, X, YFIT, YP, A) C***BEGIN PROLOGUE DP1VLU C***PURPOSE Use the coefficients generated by DPOLFT to evaluate the C polynomial fit of degree L, along with the first NDER of C its derivatives, at a specified point. C***LIBRARY SLATEC C***CATEGORY K6 C***TYPE DOUBLE PRECISION (PVALUE-S, DP1VLU-D) C***KEYWORDS CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION C***AUTHOR Shampine, L. F., (SNLA) C Davenport, S. M., (SNLA) C***DESCRIPTION C C Abstract C C The subroutine DP1VLU uses the coefficients generated by DPOLFT C to evaluate the polynomial fit of degree L, along with the first C NDER of its derivatives, at a specified point. Computationally C stable recurrence relations are used to perform this task. C C The parameters for DP1VLU are C C Input -- ALL TYPE REAL variables are DOUBLE PRECISION C L - the degree of polynomial to be evaluated. L may be C any non-negative integer which is less than or equal C to NDEG , the highest degree polynomial provided C by DPOLFT . C NDER - the number of derivatives to be evaluated. NDER C may be 0 or any positive value. If NDER is less C than 0, it will be treated as 0. C X - the argument at which the polynomial and its C derivatives are to be evaluated. C A - work and output array containing values from last C call to DPOLFT . C C Output -- ALL TYPE REAL variables are DOUBLE PRECISION C YFIT - value of the fitting polynomial of degree L at X C YP - array containing the first through NDER derivatives C of the polynomial of degree L. YP must be C dimensioned at least NDER in the calling program. C C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, C Curve fitting by polynomials in one variable, Report C SLA-74-0270, Sandia Laboratories, June 1974. C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 740601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890911 Removed unnecessary intrinsics. (WRB) C 891006 Cosmetic changes to prologue. (WRB) C 891006 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900510 Convert XERRWV calls to XERMSG calls. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DP1VLU IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER I,IC,ILO,IN,INP1,IUP,K1,K1I,K2,K3,K3P1,K3PN,K4,K4P1,K4PN, * KC,L,LM1,LP1,MAXORD,N,NDER,NDO,NDP1,NORD DOUBLE PRECISION A(*),CC,DIF,VAL,X,YFIT,YP(*) CHARACTER*8 XERN1, XERN2 C***FIRST EXECUTABLE STATEMENT DP1VLU IF (L .LT. 0) GO TO 12 NDO = MAX(NDER,0) NDO = MIN(NDO,L) MAXORD = A(1) + 0.5D0 K1 = MAXORD + 1 K2 = K1 + MAXORD K3 = K2 + MAXORD + 2 NORD = A(K3) + 0.5D0 IF (L .GT. NORD) GO TO 11 K4 = K3 + L + 1 IF (NDER .LT. 1) GO TO 2 DO 1 I = 1,NDER 1 YP(I) = 0.0D0 2 IF (L .GE. 2) GO TO 4 IF (L .EQ. 1) GO TO 3 C C L IS 0 C VAL = A(K2+1) GO TO 10 C C L IS 1 C 3 CC = A(K2+2) VAL = A(K2+1) + (X-A(2))*CC IF (NDER .GE. 1) YP(1) = CC GO TO 10 C C L IS GREATER THAN 1 C 4 NDP1 = NDO + 1 K3P1 = K3 + 1 K4P1 = K4 + 1 LP1 = L + 1 LM1 = L - 1 ILO = K3 + 3 IUP = K4 + NDP1 DO 5 I = ILO,IUP 5 A(I) = 0.0D0 DIF = X - A(LP1) KC = K2 + LP1 A(K4P1) = A(KC) A(K3P1) = A(KC-1) + DIF*A(K4P1) A(K3+2) = A(K4P1) C C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES C DO 9 I = 1,LM1 IN = L - I INP1 = IN + 1 K1I = K1 + INP1 IC = K2 + IN DIF = X - A(INP1) VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1) IF (NDO .LE. 0) GO TO 8 DO 6 N = 1,NDO K3PN = K3P1 + N K4PN = K4P1 + N 6 YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN) C C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS C DO 7 N = 1,NDO K3PN = K3P1 + N K4PN = K4P1 + N A(K4PN) = A(K3PN) 7 A(K3PN) = YP(N) 8 A(K4P1) = A(K3P1) 9 A(K3P1) = VAL C C NORMAL RETURN OR ABORT DUE TO ERROR C 10 YFIT = VAL RETURN C 11 WRITE (XERN1, '(I8)') L WRITE (XERN2, '(I8)') NORD CALL XERMSG ('SLATEC', 'DP1VLU', * 'THE ORDER OF POLYNOMIAL EVALUATION, L = ' // XERN1 // * 'REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 // * ', COMPUTED BY DPOLFT -- EXECUTION TERMINATED.', 8, 2) RETURN C 12 CALL XERMSG ('SLATEC', 'DP1VLU', + 'INVALID INPUT PARAMETER. ORDER OF POLYNOMIAL EVALUATION '// + 'REQUESTED IS NEGATIVE.', 2, 2) RETURN END SUBROUTINE DPCOEF (L, C, TC, A) C***BEGIN PROLOGUE DPCOEF C***PURPOSE Convert the DPOLFT coefficients to Taylor series form. C***LIBRARY SLATEC C***CATEGORY K1A1A2 C***TYPE DOUBLE PRECISION (PCOEF-S, DPCOEF-D) C***KEYWORDS CURVE FITTING,DATA FITTING,LEAST SQUARES,POLYNOMIAL FIT C***AUTHOR Shampine, L. F., (SNLA) C Davenport, S. M., (SNLA) C***DESCRIPTION C C Abstract C C DPOLFT computes the least squares polynomial fit of degree L as C a sum of orthogonal polynomials. DPCOEF changes this fit to its C Taylor expansion about any point C , i.e. writes the polynomial C as a sum of powers of (X-C). Taking C=0. gives the polynomial C in powers of X, but a suitable non-zero C often leads to C polynomials which are better scaled and more accurately evaluated. C C The parameters for DPCOEF are C C INPUT -- All TYPE REAL variables are DOUBLE PRECISION C L - Indicates the degree of polynomial to be changed to C its Taylor expansion. To obtain the Taylor C coefficients in reverse order, input L as the C negative of the degree desired. The absolute value C of L must be less than or equal to NDEG, the highest C degree polynomial fitted by DPOLFT . C C - The point about which the Taylor expansion is to be C made. C A - Work and output array containing values from last C call to DPOLFT . C C OUTPUT -- All TYPE REAL variables are DOUBLE PRECISION C TC - Vector containing the first LL+1 Taylor coefficients C where LL=ABS(L). If L.GT.0 , the coefficients are C in the usual Taylor series order, i.e. C P(X) = TC(1) + TC(2)*(X-C) + ... + TC(N+1)*(X-C)**N C If L .LT. 0, the coefficients are in reverse order, C i.e. C P(X) = TC(1)*(X-C)**N + ... + TC(N)*(X-C) + TC(N+1) C C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, C Curve fitting by polynomials in one variable, Report C SLA-74-0270, Sandia Laboratories, June 1974. C***ROUTINES CALLED DP1VLU C***REVISION HISTORY (YYMMDD) C 740601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891006 Cosmetic changes to prologue. (WRB) C 891006 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DPCOEF C INTEGER I,L,LL,LLP1,LLP2,NEW,NR DOUBLE PRECISION A(*),C,FAC,SAVE,TC(*) C***FIRST EXECUTABLE STATEMENT DPCOEF LL = ABS(L) LLP1 = LL + 1 CALL DP1VLU (LL,LL,C,TC(1),TC(2),A) IF (LL .LT. 2) GO TO 2 FAC = 1.0D0 DO 1 I = 3,LLP1 FAC = FAC*(I-1) 1 TC(I) = TC(I)/FAC 2 IF (L .GE. 0) GO TO 4 NR = LLP1/2 LLP2 = LL + 2 DO 3 I = 1,NR SAVE = TC(I) NEW = LLP2 - I TC(I) = TC(NEW) 3 TC(NEW) = SAVE 4 RETURN END SUBROUTINE DPOLFT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A) C***BEGIN PROLOGUE DPOLFT C***PURPOSE Fit discrete data in a least squares sense by polynomials C in one variable. C***LIBRARY SLATEC C***CATEGORY K1A1A2 C***TYPE DOUBLE PRECISION (POLFIT-S, DPOLFT-D) C***KEYWORDS CURVE FITTING, DATA FITTING,LEAST SQUARES,POLYNOMIAL FIT C***AUTHOR Shampine, L. F., (SNLA) C Davenport, S. M., (SNLA) C Huddleston, R. E., (SNLL) C***DESCRIPTION C C Abstract C C Given a collection of points X(I) and a set of values Y(I) which C correspond to some function or measurement at each of the X(I), C subroutine DPOLFT computes the weighted least-squares polynomial C fits of all degrees up to some degree either specified by the user C or determined by the routine. The fits thus obtained are in C orthogonal polynomial form. Subroutine DP1VLU may then be C called to evaluate the fitted polynomials and any of their C derivatives at any point. The subroutine DPCOEF may be used to C express the polynomial fits as powers of (X-C) for any specified C point C. C C The parameters for DPOLFT are C C Input -- All TYPE REAL variables are DOUBLE PRECISION C N - the number of data points. The arrays X, Y and W C must be dimensioned at least N (N .GE. 1). C X - array of values of the independent variable. These C values may appear in any order and need not all be C distinct. C Y - array of corresponding function values. C W - array of positive values to be used as weights. If C W(1) is negative, DPOLFT will set all the weights C to 1.0, which means unweighted least squares error C will be minimized. To minimize relative error, the C user should set the weights to: W(I) = 1.0/Y(I)**2, C I = 1,...,N . C MAXDEG - maximum degree to be allowed for polynomial fit. C MAXDEG may be any non-negative integer less than N. C Note -- MAXDEG cannot be equal to N-1 when a C statistical test is to be used for degree selection, C i.e., when input value of EPS is negative. C EPS - specifies the criterion to be used in determining C the degree of fit to be computed. C (1) If EPS is input negative, DPOLFT chooses the C degree based on a statistical F test of C significance. One of three possible C significance levels will be used: .01, .05 or C .10. If EPS=-1.0 , the routine will C automatically select one of these levels based C on the number of data points and the maximum C degree to be considered. If EPS is input as C -.01, -.05, or -.10, a significance level of C .01, .05, or .10, respectively, will be used. C (2) If EPS is set to 0., DPOLFT computes the C polynomials of degrees 0 through MAXDEG . C (3) If EPS is input positive, EPS is the RMS C error tolerance which must be satisfied by the C fitted polynomial. DPOLFT will increase the C degree of fit until this criterion is met or C until the maximum degree is reached. C C Output -- All TYPE REAL variables are DOUBLE PRECISION C NDEG - degree of the highest degree fit computed. C EPS - RMS error of the polynomial of degree NDEG . C R - vector of dimension at least NDEG containing values C of the fit of degree NDEG at each of the X(I) . C Except when the statistical test is used, these C values are more accurate than results from subroutine C DP1VLU normally are. C IERR - error flag with the following possible values. C 1 -- indicates normal execution, i.e., either C (1) the input value of EPS was negative, and the C computed polynomial fit of degree NDEG C satisfies the specified F test, or C (2) the input value of EPS was 0., and the fits of C all degrees up to MAXDEG are complete, or C (3) the input value of EPS was positive, and the C polynomial of degree NDEG satisfies the RMS C error requirement. C 2 -- invalid input parameter. At least one of the input C parameters has an illegal value and must be corrected C before DPOLFT can proceed. Valid input results C when the following restrictions are observed C N .GE. 1 C 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0. C 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0. C W(1)=-1.0 or W(I) .GT. 0., I=1,...,N . C 3 -- cannot satisfy the RMS error requirement with a C polynomial of degree no greater than MAXDEG . Best C fit found is of degree MAXDEG . C 4 -- cannot satisfy the test for significance using C current value of MAXDEG . Statistically, the C best fit found is of order NORD . (In this case, C NDEG will have one of the values: MAXDEG-2, C MAXDEG-1, or MAXDEG). Using a higher value of C MAXDEG may result in passing the test. C A - work and output array having at least 3N+3MAXDEG+3 C locations C C Note - DPOLFT calculates all fits of degrees up to and including C NDEG . Any or all of these fits can be evaluated or C expressed as powers of (X-C) using DP1VLU and DPCOEF C after just one call to DPOLFT . C C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, C Curve fitting by polynomials in one variable, Report C SLA-74-0270, Sandia Laboratories, June 1974. C***ROUTINES CALLED DP1VLU, XERMSG C***REVISION HISTORY (YYMMDD) C 740601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 891006 Cosmetic changes to prologue. (WRB) C 891006 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900911 Added variable YP to DOUBLE PRECISION declaration. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C 920527 Corrected erroneous statements in DESCRIPTION. (WRB) C***END PROLOGUE DPOLFT INTEGER I,IDEGF,IERR,J,JP1,JPAS,K1,K1PJ,K2,K2PJ,K3,K3PI,K4, * K4PI,K5,K5PI,KSIG,M,MAXDEG,MOP1,NDEG,NDER,NFAIL DOUBLE PRECISION TEMD1,TEMD2 DOUBLE PRECISION A(*),DEGF,DEN,EPS,ETST,F,FCRIT,R(*),SIG,SIGJ, * SIGJM1,SIGPAS,TEMP,X(*),XM,Y(*),YP,W(*),W1,W11 DOUBLE PRECISION CO(4,3) SAVE CO DATA CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2), 1 CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3), 2 CO(4,3)/-13.086850D0,-2.4648165D0,-3.3846535D0,-1.2973162D0, 3 -3.3381146D0,-1.7812271D0,-3.2578406D0,-1.6589279D0, 4 -1.6282703D0,-1.3152745D0,-3.2640179D0,-1.9829776D0/ C***FIRST EXECUTABLE STATEMENT DPOLFT M = ABS(N) IF (M .EQ. 0) GO TO 30 IF (MAXDEG .LT. 0) GO TO 30 A(1) = MAXDEG MOP1 = MAXDEG + 1 IF (M .LT. MOP1) GO TO 30 IF (EPS .LT. 0.0D0 .AND. M .EQ. MOP1) GO TO 30 XM = M ETST = EPS*EPS*XM IF (W(1) .LT. 0.0D0) GO TO 2 DO 1 I = 1,M IF (W(I) .LE. 0.0D0) GO TO 30 1 CONTINUE GO TO 4 2 DO 3 I = 1,M 3 W(I) = 1.0D0 4 IF (EPS .GE. 0.0D0) GO TO 8 C C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST C FOR CHOOSING DEGREE OF POLYNOMIAL FIT C IF (EPS .GT. (-.55D0)) GO TO 5 IDEGF = M - MAXDEG - 1 KSIG = 1 IF (IDEGF .LT. 10) KSIG = 2 IF (IDEGF .LT. 5) KSIG = 3 GO TO 8 5 KSIG = 1 IF (EPS .LT. (-.03D0)) KSIG = 2 IF (EPS .LT. (-.07D0)) KSIG = 3 C C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING C 8 K1 = MAXDEG + 1 K2 = K1 + MAXDEG K3 = K2 + MAXDEG + 2 K4 = K3 + M K5 = K4 + M DO 9 I = 2,K4 9 A(I) = 0.0D0 W11 = 0.0D0 IF (N .LT. 0) GO TO 11 C C UNCONSTRAINED CASE C DO 10 I = 1,M K4PI = K4 + I A(K4PI) = 1.0D0 10 W11 = W11 + W(I) GO TO 13 C C CONSTRAINED CASE C 11 DO 12 I = 1,M K4PI = K4 + I 12 W11 = W11 + W(I)*A(K4PI)**2 C C COMPUTE FIT OF DEGREE ZERO C 13 TEMD1 = 0.0D0 DO 14 I = 1,M K4PI = K4 + I TEMD1 = TEMD1 + W(I)*Y(I)*A(K4PI) 14 CONTINUE TEMD1 = TEMD1/W11 A(K2+1) = TEMD1 SIGJ = 0.0D0 DO 15 I = 1,M K4PI = K4 + I K5PI = K5 + I TEMD2 = TEMD1*A(K4PI) R(I) = TEMD2 A(K5PI) = TEMD2 - R(I) 15 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 J = 0 C C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERIA C IF (EPS) 24,26,27 C C INCREMENT DEGREE C 16 J = J + 1 JP1 = J + 1 K1PJ = K1 + J K2PJ = K2 + J SIGJM1 = SIGJ C C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1 C IF (J .GT. 1) A(K1PJ) = W11/W1 C C COMPUTE NEW A COEFFICIENT C TEMD1 = 0.0D0 DO 18 I = 1,M K4PI = K4 + I TEMD2 = A(K4PI) TEMD1 = TEMD1 + X(I)*W(I)*TEMD2*TEMD2 18 CONTINUE A(JP1) = TEMD1/W11 C C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS C W1 = W11 W11 = 0.0D0 DO 19 I = 1,M K3PI = K3 + I K4PI = K4 + I TEMP = A(K3PI) A(K3PI) = A(K4PI) A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP 19 W11 = W11 + W(I)*A(K4PI)**2 C C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE C PRECISION C TEMD1 = 0.0D0 DO 20 I = 1,M K4PI = K4 + I K5PI = K5 + I TEMD2 = W(I)*((Y(I)-R(I))-A(K5PI))*A(K4PI) 20 TEMD1 = TEMD1 + TEMD2 TEMD1 = TEMD1/W11 A(K2PJ+1) = TEMD1 C C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND C ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE C COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT, C THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST C SIGNIFICANT BITS ARE IN A(K5PI) . C SIGJ = 0.0D0 DO 21 I = 1,M K4PI = K4 + I K5PI = K5 + I TEMD2 = R(I) + A(K5PI) + TEMD1*A(K4PI) R(I) = TEMD2 A(K5PI) = TEMD2 - R(I) 21 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 C C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE C MAXDEG HAS BEEN REACHED C IF (EPS) 23,26,27 C C COMPUTE F STATISTICS (INPUT EPS .LT. 0.) C 23 IF (SIGJ .EQ. 0.0D0) GO TO 29 DEGF = M - J - 1 DEN = (CO(4,KSIG)*DEGF + 1.0D0)*DEGF FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN FCRIT = FCRIT*FCRIT F = (SIGJM1 - SIGJ)*DEGF/SIGJ IF (F .LT. FCRIT) GO TO 25 C C POLYNOMIAL OF DEGREE J SATISFIES F TEST C 24 SIGPAS = SIGJ JPAS = J NFAIL = 0 IF (MAXDEG .EQ. J) GO TO 32 GO TO 16 C C POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND. C 25 NFAIL = NFAIL + 1 IF (NFAIL .GE. 3) GO TO 29 IF (MAXDEG .EQ. J) GO TO 32 GO TO 16 C C RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT C EPS = 0.) C 26 IF (MAXDEG .EQ. J) GO TO 28 GO TO 16 C C SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.) C 27 IF (SIGJ .LE. ETST) GO TO 28 IF (MAXDEG .EQ. J) GO TO 31 GO TO 16 C C RETURNS C 28 IERR = 1 NDEG = J SIG = SIGJ GO TO 33 29 IERR = 1 NDEG = JPAS SIG = SIGPAS GO TO 33 30 IERR = 2 CALL XERMSG ('SLATEC', 'DPOLFT', 'INVALID INPUT PARAMETER.', 2, + 1) GO TO 37 31 IERR = 3 NDEG = MAXDEG SIG = SIGJ GO TO 33 32 IERR = 4 NDEG = JPAS SIG = SIGPAS C 33 A(K3) = NDEG C C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT C ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES C IF(EPS .GE. 0.0 .OR. NDEG .EQ. MAXDEG) GO TO 36 NDER = 0 DO 35 I = 1,M CALL DP1VLU (NDEG,NDER,X(I),R(I),YP,A) 35 CONTINUE 36 EPS = SQRT(SIG/XM) 37 RETURN END #endif /* DOUBLE_PRECISION */ C #ifdef SINGLE_PRECISION SUBROUTINE PVALUE (L, NDER, X, YFIT, YP, A) C***BEGIN PROLOGUE PVALUE C***PURPOSE Use the coefficients generated by POLFIT to evaluate the C polynomial fit of degree L, along with the first NDER of C its derivatives, at a specified point. C***LIBRARY SLATEC C***CATEGORY K6 C***TYPE SINGLE PRECISION (PVALUE-S, DP1VLU-D) C***KEYWORDS CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION C***AUTHOR Shampine, L. F., (SNLA) C Davenport, S. M., (SNLA) C***DESCRIPTION C C Written by L. F. Shampine and S. M. Davenport. C C Abstract C C The subroutine PVALUE uses the coefficients generated by POLFIT C to evaluate the polynomial fit of degree L, along with the first C NDER of its derivatives, at a specified point. Computationally C stable recurrence relations are used to perform this task. C C The parameters for PVALUE are C C Input -- C L - the degree of polynomial to be evaluated. L may be C any non-negative integer which is less than or equal C to NDEG , the highest degree polynomial provided C by POLFIT . C NDER - the number of derivatives to be evaluated. NDER C may be 0 or any positive value. If NDER is less C than 0, it will be treated as 0. C X - the argument at which the polynomial and its C derivatives are to be evaluated. C A - work and output array containing values from last C call to POLFIT . C C Output -- C YFIT - value of the fitting polynomial of degree L at X C YP - array containing the first through NDER derivatives C of the polynomial of degree L . YP must be C dimensioned at least NDER in the calling program. C C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, C Curve fitting by polynomials in one variable, Report C SLA-74-0270, Sandia Laboratories, June 1974. C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 740601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900510 Convert XERRWV calls to XERMSG calls. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE PVALUE DIMENSION YP(*),A(*) CHARACTER*8 XERN1, XERN2 C***FIRST EXECUTABLE STATEMENT PVALUE IF (L .LT. 0) GO TO 12 NDO = MAX(NDER,0) NDO = MIN(NDO,L) MAXORD = A(1) + 0.5 K1 = MAXORD + 1 K2 = K1 + MAXORD K3 = K2 + MAXORD + 2 NORD = A(K3) + 0.5 IF (L .GT. NORD) GO TO 11 K4 = K3 + L + 1 IF (NDER .LT. 1) GO TO 2 DO 1 I = 1,NDER 1 YP(I) = 0.0 2 IF (L .GE. 2) GO TO 4 IF (L .EQ. 1) GO TO 3 C C L IS 0 C VAL = A(K2+1) GO TO 10 C C L IS 1 C 3 CC = A(K2+2) VAL = A(K2+1) + (X-A(2))*CC IF (NDER .GE. 1) YP(1) = CC GO TO 10 C C L IS GREATER THAN 1 C 4 NDP1 = NDO + 1 K3P1 = K3 + 1 K4P1 = K4 + 1 LP1 = L + 1 LM1 = L - 1 ILO = K3 + 3 IUP = K4 + NDP1 DO 5 I = ILO,IUP 5 A(I) = 0.0 DIF = X - A(LP1) KC = K2 + LP1 A(K4P1) = A(KC) A(K3P1) = A(KC-1) + DIF*A(K4P1) A(K3+2) = A(K4P1) C C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES C DO 9 I = 1,LM1 IN = L - I INP1 = IN + 1 K1I = K1 + INP1 IC = K2 + IN DIF = X - A(INP1) VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1) IF (NDO .LE. 0) GO TO 8 DO 6 N = 1,NDO K3PN = K3P1 + N K4PN = K4P1 + N 6 YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN) C C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS C DO 7 N = 1,NDO K3PN = K3P1 + N K4PN = K4P1 + N A(K4PN) = A(K3PN) 7 A(K3PN) = YP(N) 8 A(K4P1) = A(K3P1) 9 A(K3P1) = VAL C C NORMAL RETURN OR ABORT DUE TO ERROR C 10 YFIT = VAL RETURN C 11 WRITE (XERN1, '(I8)') L WRITE (XERN2, '(I8)') NORD CALL XERMSG ('SLATEC', 'PVALUE', * 'THE ORDER OF POLYNOMIAL EVALUATION, L = ' // XERN1 // * ' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 // * ', COMPUTED BY POLFIT -- EXECUTION TERMINATED.', 8, 2) RETURN C 12 CALL XERMSG ('SLATEC', 'PVALUE', + 'INVALID INPUT PARAMETER. ORDER OF POLYNOMIAL EVALUATION ' // + 'REQUESTED IS NEGATIVE -- EXECUTION TERMINATED.', 2, 2) RETURN END SUBROUTINE PCOEF (L, C, TC, A) C***BEGIN PROLOGUE PCOEF C***PURPOSE Convert the POLFIT coefficients to Taylor series form. C***LIBRARY SLATEC C***CATEGORY K1A1A2 C***TYPE SINGLE PRECISION (PCOEF-S, DPCOEF-D) C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT C***AUTHOR Shampine, L. F., (SNLA) C Davenport, S. M., (SNLA) C***DESCRIPTION C C Written BY L. F. Shampine and S. M. Davenport. C C Abstract C C POLFIT computes the least squares polynomial fit of degree L as C a sum of orthogonal polynomials. PCOEF changes this fit to its C Taylor expansion about any point C , i.e. writes the polynomial C as a sum of powers of (X-C). Taking C=0. gives the polynomial C in powers of X, but a suitable non-zero C often leads to C polynomials which are better scaled and more accurately evaluated. C C The parameters for PCOEF are C C INPUT -- C L - Indicates the degree of polynomial to be changed to C its Taylor expansion. To obtain the Taylor C coefficients in reverse order, input L as the C negative of the degree desired. The absolute value C of L must be less than or equal to NDEG, the highest C degree polynomial fitted by POLFIT . C C - The point about which the Taylor expansion is to be C made. C A - Work and output array containing values from last C call to POLFIT . C C OUTPUT -- C TC - Vector containing the first LL+1 Taylor coefficients C where LL=ABS(L). If L.GT.0 , the coefficients are C in the usual Taylor series order, i.e. C P(X) = TC(1) + TC(2)*(X-C) + ... + TC(N+1)*(X-C)**N C If L .LT. 0, the coefficients are in reverse order, C i.e. C P(X) = TC(1)*(X-C)**N + ... + TC(N)*(X-C) + TC(N+1) C C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, C Curve fitting by polynomials in one variable, Report C SLA-74-0270, Sandia Laboratories, June 1974. C***ROUTINES CALLED PVALUE C***REVISION HISTORY (YYMMDD) C 740601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE PCOEF C DIMENSION A(*), TC(*) C***FIRST EXECUTABLE STATEMENT PCOEF LL = ABS(L) LLP1 = LL + 1 CALL PVALUE (LL,LL,C,TC(1),TC(2),A) IF (LL .LT. 2) GO TO 2 FAC = 1.0 DO 1 I = 3,LLP1 FAC = FAC*(I-1) 1 TC(I) = TC(I)/FAC 2 IF (L .GE. 0) GO TO 4 NR = LLP1/2 LLP2 = LL + 2 DO 3 I = 1,NR SAVE = TC(I) NEW = LLP2 - I TC(I) = TC(NEW) 3 TC(NEW) = SAVE 4 RETURN END SUBROUTINE POLFIT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A) C***BEGIN PROLOGUE POLFIT C***PURPOSE Fit discrete data in a least squares sense by polynomials C in one variable. C***LIBRARY SLATEC C***CATEGORY K1A1A2 C***TYPE SINGLE PRECISION (POLFIT-S, DPOLFT-D) C***KEYWORDS CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT C***AUTHOR Shampine, L. F., (SNLA) C Davenport, S. M., (SNLA) C Huddleston, R. E., (SNLL) C***DESCRIPTION C C Abstract C C Given a collection of points X(I) and a set of values Y(I) which C correspond to some function or measurement at each of the X(I), C subroutine POLFIT computes the weighted least-squares polynomial C fits of all degrees up to some degree either specified by the user C or determined by the routine. The fits thus obtained are in C orthogonal polynomial form. Subroutine PVALUE may then be C called to evaluate the fitted polynomials and any of their C derivatives at any point. The subroutine PCOEF may be used to C express the polynomial fits as powers of (X-C) for any specified C point C. C C The parameters for POLFIT are C C Input -- C N - the number of data points. The arrays X, Y and W C must be dimensioned at least N (N .GE. 1). C X - array of values of the independent variable. These C values may appear in any order and need not all be C distinct. C Y - array of corresponding function values. C W - array of positive values to be used as weights. If C W(1) is negative, POLFIT will set all the weights C to 1.0, which means unweighted least squares error C will be minimized. To minimize relative error, the C user should set the weights to: W(I) = 1.0/Y(I)**2, C I = 1,...,N . C MAXDEG - maximum degree to be allowed for polynomial fit. C MAXDEG may be any non-negative integer less than N. C Note -- MAXDEG cannot be equal to N-1 when a C statistical test is to be used for degree selection, C i.e., when input value of EPS is negative. C EPS - specifies the criterion to be used in determining C the degree of fit to be computed. C (1) If EPS is input negative, POLFIT chooses the C degree based on a statistical F test of C significance. One of three possible C significance levels will be used: .01, .05 or C .10. If EPS=-1.0 , the routine will C automatically select one of these levels based C on the number of data points and the maximum C degree to be considered. If EPS is input as C -.01, -.05, or -.10, a significance level of C .01, .05, or .10, respectively, will be used. C (2) If EPS is set to 0., POLFIT computes the C polynomials of degrees 0 through MAXDEG . C (3) If EPS is input positive, EPS is the RMS C error tolerance which must be satisfied by the C fitted polynomial. POLFIT will increase the C degree of fit until this criterion is met or C until the maximum degree is reached. C C Output -- C NDEG - degree of the highest degree fit computed. C EPS - RMS error of the polynomial of degree NDEG . C R - vector of dimension at least NDEG containing values C of the fit of degree NDEG at each of the X(I) . C Except when the statistical test is used, these C values are more accurate than results from subroutine C PVALUE normally are. C IERR - error flag with the following possible values. C 1 -- indicates normal execution, i.e., either C (1) the input value of EPS was negative, and the C computed polynomial fit of degree NDEG C satisfies the specified F test, or C (2) the input value of EPS was 0., and the fits of C all degrees up to MAXDEG are complete, or C (3) the input value of EPS was positive, and the C polynomial of degree NDEG satisfies the RMS C error requirement. C 2 -- invalid input parameter. At least one of the input C parameters has an illegal value and must be corrected C before POLFIT can proceed. Valid input results C when the following restrictions are observed C N .GE. 1 C 0 .LE. MAXDEG .LE. N-1 for EPS .GE. 0. C 0 .LE. MAXDEG .LE. N-2 for EPS .LT. 0. C W(1)=-1.0 or W(I) .GT. 0., I=1,...,N . C 3 -- cannot satisfy the RMS error requirement with a C polynomial of degree no greater than MAXDEG. Best C fit found is of degree MAXDEG . C 4 -- cannot satisfy the test for significance using C current value of MAXDEG . Statistically, the C best fit found is of order NORD . (In this case, C NDEG will have one of the values: MAXDEG-2, C MAXDEG-1, or MAXDEG). Using a higher value of C MAXDEG may result in passing the test. C A - work and output array having at least 3N+3MAXDEG+3 C locations C C Note - POLFIT calculates all fits of degrees up to and including C NDEG . Any or all of these fits can be evaluated or C expressed as powers of (X-C) using PVALUE and PCOEF C after just one call to POLFIT . C C***REFERENCES L. F. Shampine, S. M. Davenport and R. E. Huddleston, C Curve fitting by polynomials in one variable, Report C SLA-74-0270, Sandia Laboratories, June 1974. C***ROUTINES CALLED PVALUE, XERMSG C***REVISION HISTORY (YYMMDD) C 740601 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890531 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 920501 Reformatted the REFERENCES section. (WRB) C 920527 Corrected erroneous statements in DESCRIPTION. (WRB) C***END PROLOGUE POLFIT DOUBLE PRECISION TEMD1,TEMD2 DIMENSION X(*), Y(*), W(*), R(*), A(*) DIMENSION CO(4,3) SAVE CO DATA CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2), 1 CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3), 2 CO(4,3)/-13.086850,-2.4648165,-3.3846535,-1.2973162, 3 -3.3381146,-1.7812271,-3.2578406,-1.6589279, 4 -1.6282703,-1.3152745,-3.2640179,-1.9829776/ C***FIRST EXECUTABLE STATEMENT POLFIT M = ABS(N) IF (M .EQ. 0) GO TO 30 IF (MAXDEG .LT. 0) GO TO 30 A(1) = MAXDEG MOP1 = MAXDEG + 1 IF (M .LT. MOP1) GO TO 30 IF (EPS .LT. 0.0 .AND. M .EQ. MOP1) GO TO 30 XM = M ETST = EPS*EPS*XM IF (W(1) .LT. 0.0) GO TO 2 DO 1 I = 1,M IF (W(I) .LE. 0.0) GO TO 30 1 CONTINUE GO TO 4 2 DO 3 I = 1,M 3 W(I) = 1.0 4 IF (EPS .GE. 0.0) GO TO 8 C C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST C FOR CHOOSING DEGREE OF POLYNOMIAL FIT C IF (EPS .GT. (-.55)) GO TO 5 IDEGF = M - MAXDEG - 1 KSIG = 1 IF (IDEGF .LT. 10) KSIG = 2 IF (IDEGF .LT. 5) KSIG = 3 GO TO 8 5 KSIG = 1 IF (EPS .LT. (-.03)) KSIG = 2 IF (EPS .LT. (-.07)) KSIG = 3 C C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING C 8 K1 = MAXDEG + 1 K2 = K1 + MAXDEG K3 = K2 + MAXDEG + 2 K4 = K3 + M K5 = K4 + M DO 9 I = 2,K4 9 A(I) = 0.0 W11 = 0.0 IF (N .LT. 0) GO TO 11 C C UNCONSTRAINED CASE C DO 10 I = 1,M K4PI = K4 + I A(K4PI) = 1.0 10 W11 = W11 + W(I) GO TO 13 C C CONSTRAINED CASE C 11 DO 12 I = 1,M K4PI = K4 + I 12 W11 = W11 + W(I)*A(K4PI)**2 C C COMPUTE FIT OF DEGREE ZERO C 13 TEMD1 = 0.0D0 DO 14 I = 1,M K4PI = K4 + I TEMD1 = TEMD1 + DBLE(W(I))*DBLE(Y(I))*DBLE(A(K4PI)) 14 CONTINUE TEMD1 = TEMD1/DBLE(W11) A(K2+1) = TEMD1 SIGJ = 0.0 DO 15 I = 1,M K4PI = K4 + I K5PI = K5 + I TEMD2 = TEMD1*DBLE(A(K4PI)) R(I) = TEMD2 A(K5PI) = TEMD2 - DBLE(R(I)) 15 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 J = 0 C C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES DEGREE SELECTION CRITERION C IF (EPS) 24,26,27 C C INCREMENT DEGREE C 16 J = J + 1 JP1 = J + 1 K1PJ = K1 + J K2PJ = K2 + J SIGJM1 = SIGJ C C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1 C IF (J .GT. 1) A(K1PJ) = W11/W1 C C COMPUTE NEW A COEFFICIENT C TEMD1 = 0.0D0 DO 18 I = 1,M K4PI = K4 + I TEMD2 = A(K4PI) TEMD1 = TEMD1 + DBLE(X(I))*DBLE(W(I))*TEMD2*TEMD2 18 CONTINUE A(JP1) = TEMD1/DBLE(W11) C C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS C W1 = W11 W11 = 0.0 DO 19 I = 1,M K3PI = K3 + I K4PI = K4 + I TEMP = A(K3PI) A(K3PI) = A(K4PI) A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP 19 W11 = W11 + W(I)*A(K4PI)**2 C C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE C PRECISION C TEMD1 = 0.0D0 DO 20 I = 1,M K4PI = K4 + I K5PI = K5 + I TEMD2 = DBLE(W(I))*DBLE((Y(I)-R(I))-A(K5PI))*DBLE(A(K4PI)) 20 TEMD1 = TEMD1 + TEMD2 TEMD1 = TEMD1/DBLE(W11) A(K2PJ+1) = TEMD1 C C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND C ACCUMULATE SUM OF SQUARES OF ERRORS. THE POLYNOMIAL EVALUATIONS ARE C COMPUTED AND STORED IN EXTENDED PRECISION. FOR THE I-TH DATA POINT, C THE MOST SIGNIFICANT BITS ARE STORED IN R(I) , AND THE LEAST C SIGNIFICANT BITS ARE IN A(K5PI) . C SIGJ = 0.0 DO 21 I = 1,M K4PI = K4 + I K5PI = K5 + I TEMD2 = DBLE(R(I)) + DBLE(A(K5PI)) + TEMD1*DBLE(A(K4PI)) R(I) = TEMD2 A(K5PI) = TEMD2 - DBLE(R(I)) 21 SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2 C C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE C MAXDEG HAS BEEN REACHED C IF (EPS) 23,26,27 C C COMPUTE F STATISTICS (INPUT EPS .LT. 0.) C 23 IF (SIGJ .EQ. 0.0) GO TO 29 DEGF = M - J - 1 DEN = (CO(4,KSIG)*DEGF + 1.0)*DEGF FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN FCRIT = FCRIT*FCRIT F = (SIGJM1 - SIGJ)*DEGF/SIGJ IF (F .LT. FCRIT) GO TO 25 C C POLYNOMIAL OF DEGREE J SATISFIES F TEST C 24 SIGPAS = SIGJ JPAS = J NFAIL = 0 IF (MAXDEG .EQ. J) GO TO 32 GO TO 16 C C POLYNOMIAL OF DEGREE J FAILS F TEST. IF THERE HAVE BEEN THREE C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND. C 25 NFAIL = NFAIL + 1 IF (NFAIL .GE. 3) GO TO 29 IF (MAXDEG .EQ. J) GO TO 32 GO TO 16 C C RAISE THE DEGREE IF DEGREE MAXDEG HAS NOT YET BEEN REACHED (INPUT C EPS = 0.) C 26 IF (MAXDEG .EQ. J) GO TO 28 GO TO 16 C C SEE IF RMS ERROR CRITERION IS SATISFIED (INPUT EPS .GT. 0.) C 27 IF (SIGJ .LE. ETST) GO TO 28 IF (MAXDEG .EQ. J) GO TO 31 GO TO 16 C C RETURNS C 28 IERR = 1 NDEG = J SIG = SIGJ GO TO 33 29 IERR = 1 NDEG = JPAS SIG = SIGPAS GO TO 33 30 IERR = 2 CALL XERMSG ('SLATEC', 'POLFIT', 'INVALID INPUT PARAMETER.', 2, + 1) GO TO 37 31 IERR = 3 NDEG = MAXDEG SIG = SIGJ GO TO 33 32 IERR = 4 NDEG = JPAS SIG = SIGPAS C 33 A(K3) = NDEG C C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT C ALL THE DATA POINTS IF R DOES NOT ALREADY CONTAIN THESE VALUES C IF(EPS .GE. 0.0 .OR. NDEG .EQ. MAXDEG) GO TO 36 NDER = 0 DO 35 I = 1,M CALL PVALUE (NDEG,NDER,X(I),R(I),YP,A) 35 CONTINUE 36 EPS = SQRT(SIG/XM) 37 RETURN END #endif /* SINGLE_PRECISION */