C*********************************************************************** C SUBROUTINE TYPE12(TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C C ---------------------------------------------------------------------- C C TYPE 12 : COOLING and DEHUMIDIFYING COIL C C C "FORTRAN IV PROGRAM TO SIMULATE COOLING AND DEHUMIDIFYING C FINNED-TUBE MULTI-ROW HEAT EXCHANGERS", A.H. ELMAHDY AND C G.P. MITALAS, COMPUTER PROGRAM NO. 43, DIVISION OF BUILDING C RESEARCH, NATIONAL RESEARCH COUNCIL OF CANADA, OTTAWA, C MARCH 1977. C C ADAPTED FOR USE WITH MODSIM, JANUARY 1984, BY D.R. CLARK, C CENTER FOR BUILDING TECHNOLOGY, NBS. C C DYNAMICS ADDED TO STEADY STATE MODEL, THOUGH WET/DRY CONDITION C OF COIL IS DETERMINED FROM STEADY STATE CONDITIONS. C SUBROUTINES BESI AND BESK REWRITTEN IN FORTRAN 77 C NUMEROUS POINTLESS CHANGES MADE IN NOMENCLATURE C C----------------------------------------------------------------------- C C THIS PROGRAM MODELS CIRCULAR FINNED OR CONTINUOUS FINNED HEAT C EXCHANGERS WITH FOUR OR MORE ROWS IN COUNTERFLOW CROSSFLOW C CONFIGURATION. C C INPUTS: C C 1) FMR: WATER MASS FLOW RATE (KG/S) C 2) TWI: ENTERING WATER TEMPERATURE (C) C 3) FMA: DRY AIR MASS FLOW RATE (KG/S) C 4) TAIDB: ENTERING AIR DRY BULB TEMP. (C) C 5) WAI: ENTERING AIR HUMIDITY RATIO (-) C 6) TWO EXIT WATER TEMP. [FROM OUT(1)] (C) C 7) TAO EXIT AIR TEMP. [FROM OUT(2)] (C) C 8) WAO EXIT HUMID. RATIO [FROM OUT(3)] (-) C C OUTPUTS (ALL STEADY STATE VALUES): C C 1) TWO: LEAVING WATER TEMPERATURE (C) C 2) TAO: LEAVING AIR DRY BULB TEMP. (C) C 3) WAO: LEAVING AIR HUMIDITY RATIO (-) C 4) QT: TOTAL COOLING LOAD (KW) C 5) QS: SENSIBLE COOLING LOAD (KW) C 6) FWET: WET SURFACE AREA AS A FRACTION C OF TOTAL SURFACE AREA (-) C C PARAMETERS: C C 1) ICOIL: COIL TYPE: C 0 = FLAT CONTINUOUS FINS C 1 = CIRCULAR FINS (-) C 2) AP: PRIMARY SURFACE AREA (SQ. M) C 3) AS: SECONDARY SURFACE AREA (SQ. M) C 4) AI: INTERNAL SURFACE AREA (SQ. M) C 5) FLFA: MIN. FLOW AREA/FACE AREA (-) C 6) FCON: FIN MATERIAL CONDUCTIVITY (KW/M.K) C 7) FA: COIL FACE AREA (SQ. M) C 8) FPI: NUMBER OF FINS PER CENTIMETER (1/CM) C 9) NPR: NUMBER OF TUBES PER ROW (-) C 10) NOR: NUMBER OF ROWS (-) C 11) TDO: OUTSIDE TUBE DIAMETER (M) C 12) TDI: INSIDE TUBE DIAMETER (M) C 13) FTH: FIN THICKNESS (M) C 14) CM: THERMAL CAPACITANCE OF METAL (KJ/K) C 15) SL: TUBE LONGITUDINAL SPACING (M) C 16) FD: FIN DIAMETER (IF ICOIL = 1) (M) C 17) COID: COIL DEPTH (M) C 18) TCON: TUBE THERMAL CONDUCTIVITY (KW/M.K) C C*********************************************************************** DOUBLE PRECISION XIN,OUT REAL PAR,COF,S INTEGER INFO DIMENSION XIN(8),OUT(9),PAR(18),INFO(15),COF(5) C REAL PAR,SAVED,COF C INTEGER IOSTAT C LOGICAL WFLO,AFLO,WTEM,ATEM,AHUM,FRZDE C DIMENSION XIN(8),OUT(6),PAR(18),SAVED(12),COF(5) C DIMENSION IOSTAT(8) C* COMMON /CHRONO/ TIME,TSTEP,TTIME,TMIN,ITIME C* COMMON /SOSCOM/ RTOLX,ATOLX,XTOL COMMON /SIM/ TIME0,TFINAL,DELT,IWARN COMMON /STORE/ NSTORE,IAV,S(5000) DATA PI/3.141592653/ DATA FOULF/0.0003/ DATA SIACON,SILCON,SIKCON,SIQCON,SIMCON,SIFCON/10.76391, & 3.28084,577.789,3412.14,7936.68,2.54/ DATA SIT1CN,SIT2CN/1.8,32./,SIJCON/0.526543/ DATA VISC,CP,H0,SLOPE/.0437,.243,-7.55,.557/ DATA RHOA,RHOW/0.075,62.4/ DATA ADP0,ADP1,ADP2,ADP3/11.471,6875.333,-280479.327,5025099.04/ IF (INFO(7) .EQ. 0) THEN TWOI=OUT(7) TAOI=OUT(8) WAOI=OUT(9) ENDIF IF (INFO(7) .EQ. -1) THEN INFO(10)=5 CALL TYPECK(1,INFO,8,18,0) TWOI=XIN(6)*SIT1CN+SIT2CN TAOI=XIN(7)*SIT1CN+SIT2CN WAOI=XIN(8) ENDIF ISAV=INFO(10) ICOIL = PAR(1) AP = PAR(2)*SIACON AS = PAR(3)*SIACON AI = PAR(4)*SIACON FLFA = PAR(5) FCON = PAR(6)*SIKCON FA = PAR(7)*SIACON NPR = PAR(9) NOR = PAR(10) TDO = PAR(11)*SILCON TDI = PAR(12)*SILCON FTH = PAR(13)*SILCON CM = PAR(14)*SIJCON SL = PAR(15)*SILCON FD = PAR(16)*SILCON COID = PAR(17)*SILCON TCON = PAR(18)*SIKCON IF (ICOIL.EQ.0) FD=SQRT(4.*FD*COID/(PI*NOR*NPR)) FMR = XIN(1)*SIMCON IF (FMR.LT.1.) FMR=1. TWI = XIN(2)*SIT1CN+SIT2CN IF (TWI.LT.0.) TWI=0. FMA = XIN(3)*SIMCON IF (FMA.LT.1.) FMA=1. TAIDB = XIN(4)*SIT1CN+SIT2CN WAI = XIN(5) IF (WAI.LT.0.) WAI=0. TWODYN= XIN(6)*SIT1CN+SIT2CN TAODYN= XIN(7)*SIT1CN+SIT2CN WAODYN= XIN(8) C WFLO=(IOSTAT(1).LT.-1).OR.(IOSTAT(1).EQ.2) C WTEM=(IOSTAT(2).LT.-1).OR.(IOSTAT(2).EQ.2) C AFLO=(IOSTAT(3).LT.-1).OR.(IOSTAT(3).EQ.2) C ATEM=(IOSTAT(4).LT.-1).OR.(IOSTAT(4).EQ.2) C AHUM=(IOSTAT(5).LT.-1).OR.(IOSTAT(5).EQ.2) C FRZDE=WFLO.AND.AFLO.AND.WTEM.AND.ATEM.AND.AHUM AO=AP+AS HD=4.*FA*FLFA*COID/AO SREI=AO/AI FSAO=AS/AO ROH=TDO/FD HAE=0.240*TAIDB+WAI*(1061.+0.444*TAIDB) C* IF (ITIME.LE.1) THEN IF (INFO(7).EQ.-1) THEN C INITIALIZE C1,C2, AND DRY FIN EFFICIENCY COEFFICIENTS C ON FIRST CALL OF SIMULATION. CALL SUFED(ROH,COF) DO 10 I=1,5 C 10 SAVED(I)=COF(I) 10 S(ISAV+I-1)=COF(I) FH=0.5*(FD-TDO) FPI = PAR(8)*SIFCON FS=1./(12.*FPI) C SAVED(6)=0.159*(FTH/HD)**(-0.065)*(FTH/FH)**0.141 C SAVED(7)=-0.323*(FS/FH)**.049*(FD/SL)**0.549*(FTH/FS)**(-.028) C SAVED(10)=SLOPE C SAVED(11)=H0 C SAVED(12)=0. S(ISAV+5)=0.159*(FTH/HD)**(-0.065)*(FTH/FH)**0.141 S(ISAV+6)=-0.323*(FS/FH)**.049*(FD/SL)**0.549*(FTH/FS)**(-.028) S(ISAV+9)=SLOPE S(ISAV+10)=H0 S(ISAV+11)=0. ENDIF DO 28 I=1,5 C 28 COF(I)=SAVED(I) 28 COF(I)=S(ISAV+I-1) C C1=SAVED(6) C C2=SAVED(7) C IF (TIME.GT.SAVED(12)) THEN C SAVED(8)=0.5*(SAVED(8)+SAVED(10)) C SAVED(9)=0.5*(SAVED(9)+SAVED(11)) C SAVED(12)=TIME C1=S(ISAV+5) C2=S(ISAV+6) IF (TIME .GT. S(ISAV+11)) THEN S(ISAV+7)=0.5*(S(ISAV+7)+S(ISAV+9)) S(ISAV+8)=0.5*(S(ISAV+8)+S(ISAV+10)) S(ISAV+11)=TIME ENDIF C BEGIN COOLING COIL PROPER TAFRES=SREI*(.5*(TDO-TDI)/TCON+FOULF) DIA2=(TDI*12.)**.2 V=FMR*4./(NPR*PI*TDI*TDI*RHOW*3600.) FV=FMA/(RHOA*FA) C CALCULATE ENTERING AIR DEW POINT DP=ADP0+WAI*(ADP1+WAI*(ADP2+WAI*ADP3)) C CALCULATE MASS VELOCITY VM=(1.+WAI)*FMA/(FLFA*FA) FRAT=FMR/FMA FMACP=FMA*CP RENR=HD*VM/VISC C CALC. H2O FILM COEFFICIENT HCR1=(150./DIA2)*V**0.8 C IF ENTERING AIR DP DOES NOT EXCEED TWI BY MORE THAN 5 F, C SKIP TO DRY COIL SECTION. ICM=-1 IF (DP-TWI .LE. 5.) THEN TSO=DP+1. TSE=DP+1. C (next 5 lines serve ONLY to prevent compiler warning messages) HAL=HAE TWO=TWI UA=0. FIW=0. HCW=0. C (there.) GOTO 33 ENDIF C CALCULATE HCW (WET SURFACE COEFF.) CFACT=1.425-.00051*RENR+.263*RENR*RENR*10.**(-6.) HCW=0.3*CFACT*C1*RENR**C2*VM C CALCULATE THE WET FIN EFFICIENCY FTR=TWI+0.5 AWSTR=-0.00793+0.00031*FTR+0.0000075*(FTR-53.)*(FTR-53.) SHR=0.24*(TAIDB-FTR)/(.24*TAIDB+(1061.+.444*TAIDB)*WAI-.24*FTR- & (1061.+0.444*FTR)*AWSTR) SHR=ABS(SHR) AWD=WAI-AWSTR AWD=ABS(AWD) FAI=0.5*(FD-TDO)*SQRT(2.*HCW/(FCON*FTH)) IF (TAIDB.GT.80.) THEN EFFW=EXP(-0.41718)*SHR**(0.09471)*AWD**(0.0108)*FAI**(-0.50303) ELSE EFFW=EXP(-0.3574)*SHR**(0.16081)*AWD**(0.01995)*FAI**(-0.52951) ENDIF C CALCULATE OUTLET WATER TEMP AND AIR ENTHALPY FOR WET COIL FIW=FSAO*EFFW+(1.-FSAO) TRMW=TWI C SLOPE=SAVED(8) C H0=SAVED(9) SLOPE=S(ISAV+7) H0=S(ISAV+8) 41 HCRW=HCR1*(1.+0.011*TRMW) RW=SREI/HCRW RF=(CP/SLOPE)*(1.-FIW)/(FIW*HCW) RX=(TAFRES+RW)/(RF+CP/(SLOPE*HCW)) R=RX/SLOPE CWE1=1./FMA-SLOPE/FMR UAW=AO/(SLOPE*(SREI/HCRW+TAFRES)+CP/(FIW*HCW)) UA=UAW*SLOPE E=UAW*CWE1 IF (E.GT.20.) THEN TWO=(HAE-H0-(SLOPE-FRAT)*TWI)/FRAT ELSEIF (E.LT.-20.) THEN TWO=(HAE-H0)/SLOPE ELSE CW=EXP(E) TWO=((1.-CW)*(HAE-H0)+CW*(SLOPE-FRAT)*TWI)/(SLOPE-CW*FRAT) ENDIF TRMN=(TWI+TWO)/2. HAL=HAE-(FMR/FMA)*(TWO-TWI) TSO=(TWO+R*(HAE-H0))/(1.+SLOPE*R) TSE=(TWI+R*(HAL-H0))/(1.+SLOPE*R) IF (ABS(TRMN-TRMW) .GT. 0.01) THEN TRMW=TRMN HSE=H0+SLOPE*TSE HSL=H0+SLOPE*TSO TSE=TWB(HSE) TSO=TWB(HSL) SLOPE=(HSL-HSE)/(TSO-TSE) H0=HSE-SLOPE*TSE C SAVED(10)=SLOPE C SAVED(11)=H0 S(ISAV+9)=SLOPE S(ISAV+10)=H0 GOTO 41 ENDIF 33 IF (TSO .LE. DP) THEN C COMPLETELY WET ICM=0 QT=FMA*(HAE-HAL) HAI=HAE AW=AO C COMPLETELY WET COIL RESUMES AT STATEMENT C100 ELSE IF (DP .GT. TSE) ICM=1 C CALCULATE HCD (DRY SURFACE COEFF.) HCD=0.3*C1*RENR**C2*VM C CALCULATE THE DRY FIN EFFICIENCY FAI=0.5*(FD-TDO)*SQRT(2.*HCD/(FTH*FCON)) EFFD=0. DO 200 I=1,5 EFFD=EFFD+COF(I)*FAI**(I-1) 200 CONTINUE FID=FSAO*EFFD+(1.-FSAO) USD1=1./(FID*HCD) CDE1=(1./FMACP)-1./FMR IF (ICM .EQ. -1) THEN C COMPLETELY DRY COIL TRMD=TWI 32 HCRD=HCR1*(1.+0.011*TRMD) UA=AO/(USD1+SREI/HCRD+TAFRES) E=CDE1*UA IF (E.LT.-20.) THEN BETAD=FRAT/CP ELSEIF (E.GT.20.) THEN BETAD=1. ELSE CD=EXP(E) BETAD=(FRAT*(1.-CD))/(CP-FRAT*CD) ENDIF TAODB=TAIDB+BETAD*(TWI-TAIDB) QT=FMACP*(TAIDB-TAODB) TWO=TWI+QT/FMR C CALCULATE NEW TRM TRMN=(TWI+TWO)/2. IF (ABS(TRMN-TRMD) .GT. 0.01) THEN TRMD=TRMN GOTO 32 ENDIF HAL=HAE-QT/FMA TAOWB=TWB(HAL) TAOWB=TAOWB+0.0080*(TAODB-TAOWB) WAO=(HAL-.24*TAODB)/(1061.+0.444*TAODB) FWET=0. QS=QT C COMPLETELY DRY COIL STEADY STATE CALCULATION ENDS HERE GOTO 1000 ELSE C PARTLY WET COIL TRI=DP DELT=TSO-TSE C ESTIMATE WET AREA AW=AO*(DP-TSE)/(TSO-TSE) C ESTIMATE DRY AREA C SLOPE=SAVED(8) C H0=SAVED(9) SLOPE=S(ISAV+7) H0=S(ISAV+8) ERR=0. DERR=0. AWOLD=AW ICOUNT=0 90 AD=AO-AW TRMD=(TRI+TWO)/2. HCRD=HCR1*(1.+0.011*TRMD) UAD=AD/(USD1+SREI/HCRD) TRMW=(TWI+TRI)/2. HCRW=HCR1*(1.+0.011*TRMW) RW=SREI/HCRW RF=(CP/SLOPE)*(1.-FIW)/(FIW*HCW) RX=(TAFRES+RW)/(RF+CP/(SLOPE*HCW)) R=RX/SLOPE CWE1=1./FMA-SLOPE/FMR UAW=AW/(SLOPE*(SREI/HCRW+TAFRES)+CP/(FIW*HCW)) UA=UAD+UAW*SLOPE E=CDE1*UAD IF (E.LT.-20.) THEN BETAD=FRAT/CP ELSEIF (E.GT.20.) THEN BETAD=1. ELSE CD=EXP(E) BETAD=(FRAT*(1.-CD))/(CP-FRAT*CD) ENDIF CPBETA=CP*BETAD E=CWE1*UAW IF (E.GT.20.) THEN TRI=((SLOPE-FRAT)*TWI-(HAE-H0-CPBETA*TAIDB)) & /(CPBETA-FRAT) ELSEIF (E.LT.-20.) THEN TRI=(HAE-H0-CPBETA*TAIDB)/(SLOPE-CPBETA) ELSE CW=EXP(E) TRI=(CW*(SLOPE-FRAT)*TWI+(1.-CW)*(HAE-H0-CPBETA*TAIDB))/ & (SLOPE-CW*FRAT-(1.-CW)*CPBETA) ENDIF TAI=TAIDB+BETAD*(TRI-TAIDB) HAI=HAE-CP*(TAIDB-TAI) TSI=(TRI+R*(HAI-H0))/(1.+SLOPE*R) HAL=HAI-FRAT*(TRI-TWI) TSE=(TWI+R*(HAL-H0))/(1.+SLOPE*R) QT=FMA*(HAE-HAL) TWO=TWI+QT/FMR IF (ICOUNT.GE.1) DERR=DP-TSI-ERR ERR=DP-TSI IF (ABS(ERR) .GT. 0.01) THEN IF (ICOUNT.LT.1 .OR. ABS(DERR).LT.1.E-6) THEN AWOLD=AW AW=AW+0.7*AO*ERR/(TSO-TSE) ELSE DAW=AW-AWOLD AWOLD=AW AW=AW-ERR*DAW/DERR ENDIF ICOUNT=ICOUNT+1 IF (AW.GT.AO)AW=AO IF (AW.LE.0.)AW=0.001*AO HSE=H0+SLOPE*TSE HSI=H0+SLOPE*TSI TSE=TWB(HSE) TSI=TWB(HSI) SLOPE=(HSI-HSE)/(TSI-TSE) H0=HSE-SLOPE*TSE C SAVED(10)=SLOPE C SAVED(11)=H0 S(ISAV+9)=SLOPE S(ISAV+10)=H0 GOTO 90 ENDIF ENDIF ENDIF C PARTLY OR COMPLETELY WET COIL C100 TAOWB=TWB(HAL) C CALCULATE DRY BULB TEMPERATURE C=(HCW*AW)/(CP*FMA) C=C*FIW Y=0. IF (C.LT.20.)Y=EXP(-C) HSP=HAI-(HAI-HAL)/(1.-Y) TSP=TWB(HSP) TAODB=TSP+(TAIDB-TSP)*Y WAO=(HAL-.24*TAODB)/(1061.+0.444*TAODB) QS=0.243*FMA*(TAIDB-TAODB) FWET=AW/AO 1000 CONTINUE C** END OF STEADY STATE MODEL; TACK ON SIMPLE DYNAMICS C** STEADY STATE OUTPUT IF PAR(14) .LE. 1.0 IF (CM.LE.SIJCON) THEN C DTAODT=TAODB-SIT2CN C DTWODT=TWO-SIT2CN C DWAODT=WAO TAO=TAODB TWO=TWO C IOSTAT(1)=1 C IOSTAT(2)=1 C IOSTAT(3)=1 ELSE RATE=UA/(3600.*CM) C DTAODT=RATE*(TAODB-TAODYN) C DWAODT=RATE*(WAO-WAODYN) C DTWODT=RATE*(TWO-TWODYN) TAOAA=0 TAOBB=RATE*(TAODB-TAODYN) CALL DIFFEQ(TIME,TAOAA,TAOBB,TAOI,TAOF,TAOBAR) TAO=TAOBAR WAOAA=RATE WAOBB=-RATE/WAODYN CALL DIFFEQ(TIME,WAOAA,WAOBB,WAOI,WAOF,WAOBAR) WAO=WAOBAR TWOAA=RATE TWOBB=-RATE/TWODYN CALL DIFFEQ(TIME,TWOAA,TWOBB,TWOI,TWOF,TWOBAR) TWO=TWOBAR C IOSTAT(1)=0 C IOSTAT(2)=0 C IOSTAT(3)=0 C IF ((ABS(TWO-TWODYN).LE.RTOLX*ABS(TWO)+ATOLX) .AND. C & (ABS(TAODB-TAODYN).LE.RTOLX*ABS(TAODB)+ATOLX) .AND. C & (ABS(WAO-WAODYN).LE.RTOLX*ABS(WAO)+ATOLX).AND.FRZDE) C & THEN C IOSTAT(1)=1 C IOSTAT(2)=1 C IOSTAT(3)=1 C ENDIF ENDIF C** CONVERT UNITS FOR OUTPUT (ALL COIL CONDITIONS) C ----------------Y. ZHU------------------ C IF (ABS(DTWODT).LT.1.E-6) DTWODT=0.0 ! added 2/24/94 C IF (ABS(DTAODT).LT.1.E-6) DTAODT=0.0 C IF (ABS(DWAODT).LT.1.E-6) DWAODT=0.0 C ---------------------------------------- C OUT(1)=DTWODT/SIT1CN C OUT(2)=DTAODT/SIT1CN C OUT(3)=DWAODT OUT(1)=(TWO-SIT2CN)/SIT1CN OUT(2)=(TAO-SIT2CN)/SIT1CN OUT(3)=WAO OUT(4)=QT/SIQCON OUT(5)=QS/SIQCON OUT(6)=FWET OUT(7)=TWOF OUT(8)=TAOF OUT(9)=WAOF C IOSTAT(4)=1 C IOSTAT(5)=1 C IOSTAT(6)=1 RETURN 1 END C*********************************************************************** FUNCTION TWB(H) C*********************************************************************** C THE FOLLOWING FUNCTION DEFINES WET BULB TEMPERATURE C AS A FUNCTION OF ALOG(ENTHALPY): C IF (H.LT.11.83)THEN TWB=32. ELSE ALOGH=ALOG(H) TWB=30.9185+ALOGH*(-39.682+ALOGH*(20.5841-1.758*ALOGH)) ENDIF RETURN END C*********************************************************************** SUBROUTINE SUFED(ROH,COF) C*********************************************************************** DIMENSION X(60),Y(60),COF(5) C*** ROH = (OUTSIDE TUBE DIAM.)/(EFFECTIVE FIN DIAM.) FAI=0.02 DO 10 I=1,60 FAI=FAI+0.035 R1=FAI/(1.-ROH) R2=R1*ROH RO=2.*ROH/(FAI*(1.+ROH)) CALL BESI(R1,1,R1I1,IE1) CALL BESK(R2,1,R2K1,IE2) CALL BESI(R2,1,R2I1,IE3) CALL BESK(R1,1,R1K1,IE4) CALL BESI(R2,0,R2I0,IE5) CALL BESK(R2,0,R2K0,IE6) FED=RO*(R1I1*R2K1-R2I1*R1K1)/(R2I0*R1K1+R1I1*R2K0) X(I)=FAI Y(I)=FED 10 CONTINUE CALL POLFIT(X,Y,60,1.E-05,4,COF) RETURN END C*********************************************************************** SUBROUTINE BESI(X,N,BI,IER) C ---------------------------------------------------------------------- C C*** SUBROUTINE BESI TO CALCULATE THE MODIFIED BESSEL FUNCTION C*** OF ORDER 0 TO N C*** C*** CALL BESI(X,N,BI,IER) C*** C*** X ARGUMENT OF BESSEL FUNCTION C*** N ORDER OF BESSEL FUNCTION, GREATER THAN OR EQUAL TO ZERO C*** BI RESULTANT VALUE OF I BESSEL FUNCTION C*** IER RESULTANT ERROR CODE: C*** IER = 0 NO ERROR C*** IER = 1 N .LT. 0 C*** IER = 2 X .LT. 0 C*** IER = 3 BI .LT. 10**(-30), BI IS SET TO 0 C*** IER = 4 X .GT. N & X .GT. 90, BI IS SET TO 10**38 C*** C*********************************************************************** DATA PI,TOL/3.141592653,1.E-06/ IER=0 BI=1. IF (X .EQ. 0. .AND. N .EQ. 0) RETURN IF (N .LT. 0) THEN IER=1 RETURN ELSE IF (X .LT. 0.) THEN IER = 2 RETURN ELSE IF (X .GT. 12. .AND. X .GT. N) THEN IF (X .GT. 90.) THEN IER = 4 BI=10.**30 RETURN ENDIF FN=4*N*N XX=0.125/X TERM=1. BI=1. DO 10 K=1,30 IF (ABS(TERM).LE.ABS(TOL*BI)) THEN BI=BI*EXP(X)/SQRT(2.*PI*X) RETURN ENDIF FK=(2*K-1)**2 TERM=TERM*XX*(FK-FN)/REAL(K) BI=BI+TERM 10 CONTINUE ENDIF XX=X/2. TERM=1. IF (N.GT.0) THEN DO 20 I=1,N FI=I IF (ABS(TERM).LT.1.E-30*FI/XX) THEN IER=3 BI=0. RETURN ENDIF TERM=TERM*XX/FI 20 CONTINUE ENDIF BI=TERM XX=XX*XX DO 30 K=1,1000 IF (ABS(TERM).LE.ABS(BI*TOL)) RETURN FK=K*(N+K) TERM=TERM*(XX/FK) BI=BI+TERM 30 CONTINUE RETURN END C*********************************************************************** SUBROUTINE BESK(X,N,BK,IER) C ---------------------------------------------------------------------- C C*** SUBROUTINE BESK TO COMPUTE THE K BESSEL FUNCTION FOR A GIVEN C*** ARGUMENT AND ORDER C*** C*** CALL BESK(X,N,BK,IER) C*** C*** X THE ARGUMENT OF THE K BESSEL FUNCTION DESIRED C*** N THE ORDER OF THE K BESSEL FUNCTION DESIRED C*** BK THE RESULTANT K BESSEL FUNCTION C*** IER RESULTANT ERROR CODE: C*** IER=0 NO ERROR C*** IER=1 N IS NEGATIVE C*** IER=2 X IS ZERO OR NEGATIVE C*** IER=3 X .GT. 85, BK .LT. 10**-38; BK SET TO 0. C*** IER=4 BK .GT. 10**38; BK SET TO 10**38 C*** C*** NOTE: N MUST BE GREATER THAN OR EQUAL TO ZERO C*** C*** METHOD: C*** COMPUTES ZERO ORDER AND FIRST ORDER BESSEL FUNCTIONS USING C*** SERIES APPROXIMATIONS AND THEN COMPUTES N TH ORDER FUNCTION C*** USING RECURRENCE RELATION. C*** RECURRENCE RELATION AND POLYNOMIAL APPROXIMATION TECHNIQUE C*** AS DESCRIBED BY A.J.M. HITCHCOCK, 'POLYNOMIAL APPROXIMATIONS C*** TO BESSEL FUNCTIONS OF ORDER ZERO AND ONE AND TO RELATED C*** FUNCTIONS,' M.T.A.C., V.11, 1957, PP. 86-88, AND G.N. WATSON, C*** 'A TREATISE ON THE THEORY OF BESSEL FUNCTIONS,' CAMBRIDGE C*** UNIVERSITY PRESS, 1958, P.62 C*** C*********************************************************************** DIMENSION T(12) DATA GJMAX/ 1.E38/ BK=0. G0=0. GJ=0. IF (N.LT.0) THEN IER=1 RETURN ELSE IF (X.LE.0.) THEN IER=2 RETURN ELSE IF (X.GT.85.) THEN IER=3 BK=0. RETURN ENDIF IER=0 C Use polynomial approximation if x > 1. IF (X.GT.1.) THEN A=EXP(-X) B=1./X C=SQRT(B) T(1)=B DO 10 L=2,12 T(L)=T(L-1)*B 10 CONTINUE IF (N.NE.1) THEN C Compute K0 using polynomial approximation G0=A*(1.2533141-.1566642*T(1)+.08811128*T(2)-.09139095 & *T(3)+.1344596*T(4)-.2299850*T(5)+.3792410*T(6)-.5247277 & *T(7)+.5575368*T(8)-.4262633*T(9)+.2184518*T(10) & -.06680977*T(11)+.009189383*T(12))*C IF (N .EQ. 0) THEN BK=G0 RETURN ENDIF ENDIF C Compute K1 using polynomial approximation G1=A*(1.2533141+.4699927*T(1)-.1468583*T(2)+.1280427*T(3) & -.1736432*T(4)+.2847618*T(5)-.4594342*T(6)+.6283381*T(7) & -.6632295*T(8)+.5050239*T(9)-.2581304*T(10)+.07880001*T(11) & -.01082418*T(12))*C IF (N .EQ. 1) THEN BK=G1 RETURN ENDIF ELSE C Use series expansion if x <= 1. B=X/2. A=.5772157+ALOG(B) C=B*B IF (N.NE.1) THEN C Compute K0 using series expansion G0=-A X2J=1. FACT=1. HJ=0. DO 20 J=1,6 RJ=1./FLOAT(J) X2J=X2J*C FACT=FACT*RJ*RJ HJ=HJ+RJ G0=G0+X2J*FACT*(HJ-A) 20 CONTINUE IF (N.EQ.0) THEN BK=G0 RETURN ENDIF ENDIF C Compute K1 using series expansion X2J=B FACT=1. HJ=1. G1=1./X+X2J*(.5+A-HJ) DO 30 J=2,8 X2J=X2J*C RJ=1./FLOAT(J) FACT=FACT*RJ*RJ HJ=HJ+RJ G1=G1+X2J*FACT*(.5+(A-HJ)*FLOAT(J)) 30 CONTINUE IF (N.EQ.1) THEN BK=G1 RETURN ENDIF ENDIF C From K0 and K1 compute KN using recurrence relation DO 40 J=2,N GJ=2.*(FLOAT(J)-1.)*G1/X+G0 IF (GJ-GJMAX) 33,33,32 32 IER=4 GJ=GJMAX GOTO 34 33 G0=G1 G1=GJ 40 CONTINUE 34 BK=GJ RETURN END C*********************************************************************** SUBROUTINE POLFIT(X,Y,N,TOL,LAST,COF) C --------------------------------------------------------------------- C C*** SUBROUTINE POLFIT FITS POLYNOMIAL OF ORDER FROM 1 TO LAST TO THE C*** ORDERED PAIRS OF DATA POINTS X,Y C*** C*********************************************************************** DIMENSION X(60),Y(60), A(10,10),SUMX(10),SUMY(10),COF(10) SUMX(1)=N SUMX(2)=0. SUMX(3)=0. SUMY(1)=0. SUMY(2)=0. L=LAST+1 DO 80 I=1,L 80 COF(I)=0. DO 90 I=1,N SUMX(2)=SUMX(2)+X(I) SUMX(3)=SUMX(3)+X(I)*X(I) SUMY(1)=SUMY(1)+Y(I) SUMY(2)=SUMY(2)+X(I)*Y(I) 90 CONTINUE NORD=1 91 L=NORD+1 KK=L+1 DO 101 I=1,L DO 100 J=1,L IK=J-1+I A(I,J)=SUMX(IK) 100 CONTINUE A(I,KK)=SUMY(I) 101 CONTINUE DO 140 I=1,L A(KK,I)=-1. KKK=I+1 DO 110 J=KKK,KK A(KK,J)=0. 110 CONTINUE C=1./A(1,I) DO 120 II=2,KK DO 120 J=KKK,KK A(II,J)=A(II,J)-A(1,J)*A(II,I)*C 120 CONTINUE DO 140 II=1,L DO 140 J=KKK,KK A(II,J)=A(II+1,J) 140 CONTINUE S2=0. DO 160 J=1,N S1=A(1,KK) DO 150 I=1,NORD S1=S1+A(I+1,KK)*X(J)**I 150 CONTINUE S2=S2+(S1-Y(J))*(S1-Y(J)) 160 CONTINUE B=N-L IF (S2 .GT. 0.0001) S2=SQRT(S2/B) DO 164 I=1,L J=I-1 COF(I)=A(I,KK) 164 CONTINUE IF (NORD-LAST) 170,173,173 170 IF (S2-TOL) 173,173,171 171 NORD=NORD+1 J=2*NORD SUMX(J)=0. SUMX(J+1)=0. SUMY(NORD+1)=0. DO 172 I=1,N SUMX(J)=SUMX(J)+X(I)**(J-1) SUMX(J+1)=SUMX(J+1)+X(I)**J SUMY(NORD+1)=SUMY(NORD+1)+Y(I)*X(I)**NORD 172 CONTINUE GOTO 91 173 CONTINUE RETURN END