*CMZ : 05/07/94 12.48.22 by Unknown *-- Author : * * REAL*8 FUNCTION SAKAP(N,DT,TB) ********************************************************************** * * * CALCULATE KAPITZA CONDUCTIVITY FUNCTIONS(W/M.M) * * N=0: Constant given by the input parameter KAPITZA [W/(m^2 K)] * * N=1: UNANNEALED NIOBIUM * * N=2: ANNEALED NIOBIUM * * N=3: NUCLEATE BOILING HELIUM * * N=4: SUBCOOLED HELIUM * * N=5: HELIUM EXCHANGE * C NOTE: With TBATH=1.4, N=3 starts at T-TBATH=1.98. BPG IMPLICIT REAL*8(A-H,O-Z) REAL DUMMY REAL*8 DUMMY8 IF (DT.LT.0.)THEN IF (DT.GE.-0.000005)THEN DT=0.D0 ELSE WRITE(6,1018)DT,TB 1018 FORMAT(///' WARNING: TEMP IS LESS THAN TBATH', + /' DT = ',F6.4) SAKAP=0.D0 RETURN ENDIF ENDIF C ***** IMPORTANT: Uncommment this region to allow for the C ***** Kapitza models to be used. The function C ***** always returns a constant now. C C* Constant for debugging C IF(N.EQ.0) THEN C call get_kapitza(DUMMY) C DUMMY8 = DUMMY C SAKAP = DUMMY8*DT C RETURN C ENDIF C C IF(N.EQ.3)THEN CC BPG This is suspect! I'd like 12302.7!!!!!!!!!!!!! CC SAKAP=10230.27*DT**0.45 C SAKAP=12302.7*DT**0.45 C RETURN C ENDIF C CC All below agree with HEAT93.FLX BPG 4/14/94 C IF(N.EQ.4)THEN C SAKAP=154.9267*DT**0.333*(TB/2.5)**1.6 C RETURN C ENDIF C C TP=DT/TB C C P = .25*TP*TP*TP C P = P + 1.0 + 1.5*TP + TP*TP C IF(N.EQ.1)SAKAP=170.*P*TB**3.62 C IF(N.EQ.2)SAKAP=200.0*P*(TB**4.65) C IF(N.EQ.5)SAKAP=10.**(1.126781+1.079807*LOG10(DT)) * Constant for debugging call get_kapitza(DUMMY) DUMMY8 = DUMMY SAKAP = DUMMY8*DT RETURN END *CMZ : 11/07/94 14.03.53 by Unknown *-- Author : Cao * * REAL*8 FUNCTION COND(N,T,TC) ********************************************************************** * * * CALCULATE THERMAL CONDUCTIVITY [W/(m K)] * * N=1: RRR=20 (WAH CHANG) * * N=2: RRR=40 (KBI) * * N=3: RRR=120 (HERAEUS) * * N=4: RRR=30000 (GLANDUN) * * N=5: AN OLD THERMAL CONDUCTIVITY REPORTED BY CERN * * N=6: RRR=300 (WASSERBACH) * * N=7: RRR=500 (THEORETICAL) * * N=8: RRR=1000 (THEORETICAL) * * N=9: RRR=2000 (THEORETICAL) * * N=10: 30 RRR NIOBIUM WITH LOW PHONON CONTRIBUTION * * N=11: 30 RRR NIOBIUM WITH FLAT PHONON CONTRIBUTION * * N=12: 30 RRR NIOBIUM WITH HIGH PHONON CONTRIBUTION * * N=90: 90 RRR from HEAT programs * * N=250: 250 RRR from HEAT programs * * N=400: 400 RRR from HEAT programs * * N=403: 400 RRR from HEAT programs + Level off at 500 * * N=401: 400 RRR from HEAT programs + Bigger phonon peak * * N=420: 400 RRR from HEAT programs + 2(elect. cont) * * N=440: 400 RRR from HEAT programs + 4(elect. cont) * IMPLICIT REAL*8(A-H,O-Z) REAL DUMMY REAL*8 DUMMY8 C--------------- Executable Code starts here ----------------------- * Constant for debugging IF(N.EQ.0) THEN call get_conduct(DUMMY) DUMMY8 = DUMMY COND = DUMMY8 RETURN ENDIF * RRR = 403, from HEAT program, 08 Oct 91 + level off at 500 IF(N.EQ.403)THEN X=LOG10(T) IF(X.LE.0.40)THEN Y=-5.6007+131.09*X-972.17*X**2+3623.2*X**3 + -6710.8*X**4+4889.5*X**5 COND=10**Y ELSE IF(X.LE.1.210)THEN Y=11.759-73.392*X+192.27*X**2-233.17*X**3 + +137.65*X**4-32.139*X**5 COND=10**Y ELSE IF(X.LE.1.496)THEN P=T+4.35 COND= 3965.0-150.0*P+4.0E-2*P**3 ELSE COND= 374.5+1.94*T-6.94E-3*T**2 ENDIF RETURN ENDIF * RRR = 400, from HEAT program, 08 Oct 91 IF(N.EQ.400)THEN X=LOG10(T) IF(X.LE.0.40)THEN Y=-5.6007+131.09*X-972.17*X**2+3623.2*X**3 + -6710.8*X**4+4889.5*X**5 ELSE Y=11.759-73.392*X+192.27*X**2-233.17*X**3 + +137.65*X**4-32.139*X**5 ENDIF COND=10**Y RETURN ENDIF * RRR = 400, from HEAT program, 08 Oct 91 + Bigger phonon peak C Draw a line between the ends of the phonon peak, subtract this C line from the total function, what's left is the "phonon peak" C New COND is the old line + 15("phonon peak") C This is a test for HSP BPG IF(N.EQ.401)THEN X=LOG10(T) IF(X.LE.0.40)THEN Y=-5.6007+131.09*X-972.17*X**2+3623.2*X**3 + -6710.8*X**4+4889.5*X**5 ELSE Y=11.759-73.392*X+192.27*X**2-233.17*X**3 + +137.65*X**4-32.139*X**5 ENDIF IF(X .LE. 0.4)THEN COND=15.D0*(10**Y-(17.569+4.1006*T))+ + (17.569+4.1006*T) ELSE COND=10**Y ENDIF RETURN ENDIF * RRR = 400, from HEAT program, 08 Oct 91 + 2(electronic cont) C Multiply the high temp part by 4 and add interpolating fcn. IF(N.EQ.420)THEN X=LOG10(T) IF(X.LE.0.34242268)THEN Y=-5.6007+131.09*X-972.17*X**2+3623.2*X**3 + -6710.8*X**4+4889.5*X**5 COND=10**Y ELSE IF(X.LE.0.54406804)THEN COND=111.51-58.983*T+1.4802*T**2+4.2332*T**3 ELSE Y=11.759-73.392*X+192.27*X**2-233.17*X**3 + +137.65*X**4-32.139*X**5 COND=2.*(10**Y) ENDIF RETURN ENDIF * RRR = 400, from HEAT program, 08 Oct 91 + 4(electronic cont) C Multiply the high temp part by 4 and add interpolating fcn. IF(N.EQ.440)THEN X=LOG10(T) IF(X.LE.0.34242268)THEN Y=-5.6007+131.09*X-972.17*X**2+3623.2*X**3 + -6710.8*X**4+4889.5*X**5 COND=10**Y ELSE IF(X.LE.0.54406804)THEN COND=377.33-368.53*T+104.83*T**2-3.7576*T**3 ELSE Y=11.759-73.392*X+192.27*X**2-233.17*X**3 + +137.65*X**4-32.139*X**5 COND=4.*(10**Y) ENDIF RETURN ENDIF * RRR = 90, from HEAT program, 08 Oct 91 IF(N.EQ.90)THEN X=LOG10(T) IF(X.LE.0.44)THEN Y=-9.0695+161.05*X-1015.6*X**2+3173.5*X**3 + -4882.5*X**4+2951.7*X**5 ELSE Y=-12.074+106.50*X-334.67*X**2+509.33*X**3 + -371.17*X**4+104.12*X**5 ENDIF COND=10**Y RETURN ENDIF * RRR = 250, from HEAT program, 08 Oct 91 IF(N.EQ.250)THEN X=LOG10(T) IF(X.LE.0.48)THEN Y=-6.8734+91.059*X-276.76*X**2+41.596*X**3 + +919.88*X**4-986.54*X**5 ELSE Y=-22.482+173.46*X-511.72*X**2+756.08*X**3 + -548.63*X**4+156.23*X**5 ENDIF COND=10**Y RETURN ENDIF * RRR = 300, Wasserbach, 08 Oct 91 IF(N.EQ.6) THEN X=T*1.533 IF(X.LT.4.56) THEN X=LOG10(X) Y=6.267085-8.992381*X+64.15068*X*X-63.22761*X**3 ELSEIF(X.LT.11.10) THEN Y=107.4531-64.69516*X+16.12723*X*X- * 1.887882*X**3+0.1064409*X**4-2.3366571E-3*X**5 ELSE Y=13.25630+0.2958537*X ENDIF COND=10**(Y/6.76) RETURN ENDIF * RRR=500 OR 1000 OR 2000 IF(N.EQ.7.OR.N.EQ.8.OR.N.EQ.9)THEN IF(N.EQ.7) RRR=500 IF(N.EQ.8) RRR=1000 IF(N.EQ.9) RRR=2000 X=T/9.2*14.5 IF(X.LT.3.39) THEN COND=.033935*0.1674*RRR*2.22 RETURN ELSE X=LOG10(X) Y=-4.895598+12.75957*X-8.926826*X*X+2.08717*X**3 Y=10**Y Y=Y/13.85 COND=Y*0.1674*RRR*T RETURN ENDIF ENDIF * LO CURVE. VALID BETWEEN 1.5 TO 4 KELVIN IF(N.EQ.10) THEN IF(T.LT.1.5) THEN COND=.385*(((T/1.5)**3)) ELSE X=LOG10(T) Y=5.281229E-2+3.845072*X-1.244191*X*X COND=10.**Y/10. ENDIF RETURN ENDIF * FLAT CURVE IF(N.EQ.11) THEN IF(T.LT.1.5) THEN COND=5.6*((T/1.5)**3) ELSE X=LOG10(T) Y=-9.4661156E-2+8.300755*X-23.56781*X*X+20.54456*X**3 COND=10.**Y ENDIF RETURN ENDIF * HIGH CURVE IF(N.EQ.12) THEN X=LOG10(T) Y=0.9595242+10.22286*X-38.84547*X*X+36.29954*X**3 COND=10**Y RETURN ENDIF * CERN'S WRONG THERMAL CONDUCTIVITY IF(N.EQ.5) THEN IF(T.GT.6.) THEN COND=5.5 ELSE X=LOG10(T*1.775) Y=7.744585-56.03619*X+114.1993*X**2-55.65905*X**3 COND=10.**(Y/14.2) ENDIF RETURN ENDIF * THERMAL CONDUCTIVITY OF GLANDUN ETAL NB ID USED IF(N.EQ.4) THEN X=2.05*LOG(T) IF(X.LT..90)Y=1.8889*X+2.0 IF(X.GE..9.AND.X.LT.2.0)Y=4.3333*X-.2167 IF(X.GE.2.0.AND.X.LT.3.0)Y=-1.8*X*X+10.6*X-5.55 IF(X.GE.3.0)Y=.10*X+9.75 COND=13.0*EXP(Y/2.05) ELSE X=LOG10(T)*21. * WAN CHANG THERMAL CONDUCTIVITY IF(N.EQ.1) THEN IF(X.GT.10.4) THEN Y=-14.43890+1.618171*X+1.4486675E-2*X*X-1.641114E-3*X**3 ELSEIF(T.LT.2.) THEN Y=0. ELSE Y=-3.390097+.5364078*X ENDIF ENDIF * KBI THERMAL CONDUCTIVITY IF(N.EQ.2) THEN IF(X.GT.7.43) THEN Y=-81.503444+36.71027*X-5.908534*X*X+.4527901*X*X*X Y=Y-1.6412885E-2*X*X*X*X+2.2711556E-4*X*X*X*X*X ELSE Y=-1.390402+0.9865952*X ENDIF ENDIF * HERAEUS THERMAL CONDUCTIVITY IF(N.EQ.3) THEN IF(X.GT.8.9)THEN Y=-63.14782+29.95032*X-5.189883*X*X+0.4380147*X*X*X Y=Y-1.744542E-2*X*X*X*X+2.6300667E-4*X**5 ELSE Y=-1.26431+1.016857*X-1.7964604E-2*X**2 ENDIF ENDIF COND=10.**(Y/7.54) IF(COND.LT.1.) COND=1. ENDIF END *CMZ : 22/07/94 13.49.13 by Unknown *-- Author : * * REAL*8 FUNCTION HECA(NSP,T,TC) ********************************************************************** * * * CALCULATION OF (inverse) HEAT CAPACITY OF NIOBIUM * * CN=C*T+B*T**3; CS=A*T**3+B*T**3 * * A+B=36.2; B=8.25; C=698.0 * * Corrections to heat capacity put in 29 Oct 91; JHG * C NSP = 1 Original model supplied by Cao, et.al with THEAT code. C As fixed by JHG IMPLICIT REAL*8(A-H,O-Z) REAL DUMMY REAL*8 DUMMY8 * Constant value for debugging IF(NSP.EQ.0)THEN call get_heatcap(DUMMY) DUMMY8 = DUMMY HECA = 1.0D0/DUMMY8 ENDIF IF(NSP.EQ.1)THEN C This is the original function C Looks like this function of C has units (J/(K.m^3) C BPG Not identical but compares well to VanDerHoeven C & Keesom. IF(T.GT.TC)THEN HECA=1.D0/(698.D0*T+8.25D0*T*T*T) ELSE IF (T.GT.4.D0)THEN HECA=1.D0/(36.2D0*T*T*T) ELSE HECA=1.D0/ + (-185.03D0+310.84D0*T-231.24D0*T*T+77.128D0*T*T*T) ENDIF ENDIF END *CMZ : 05/07/94 12.48.32 by Unknown *-- Author : * * REAL*8 FUNCTION SURE(T) ********************************************************************** * * * CALCULATION OF BCS SURFACE RESISTANCE OF NIOBIUM * * See for example CLNS-79/434 p3. * T - Absolute temperature * FREQUENCY F (GHZ) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONT/HMAX,HM,RED,RN,TBATH,F,EX2,C1,RS,TC, + TCS,FAC,BCN,NSP,NT,NK,NH IF(T.GT.TC) THEN SURE=RN ELSE FN=F/2.856D0 C 9.22 is Tc0, i.e. Tc in zero magnetic field (labelled TC here) C Want TCS? No! HSP says use no H-dependence in Gap !!!!!!!! TT=T/9.22 C ***** IMPORTANT : This value has been made constant. Change C ***** this back for a real simulation. C SURE=FAC*1.61E-4*FN*FN/T*LOG(16.*T/FN)*EXP(-17.2*SQRT C * (COS(3.1415926*TT*TT/2.))/T)+RS SURE = 1.55D-6 ENDIF C DEBUG C write(*,*)' TEMP = ',T,' SURE = ',SURE END