FUNCTION RCU(RRRR,TTI,BBF)
C** function to calculate resistance of Cu as a function of 
C   temp and RRR.  Also calculates magneto resistance for BF
C   in Tesla.  Uses POLFIT to interpolate. 
      implicit double precision (a-h,o-z)
      real rrrr,tti,bbf
      dimension b(4),x(4),y(4)
      data r273,rup,r373,rupp/1.553,4.28d-3,1.981,3.8d-3/
C  *** points where curves merge:
C   RRR1000=500 at 50
C   =200 at 60
C   =100 at 80
C   =50 at 150
C   =20 at 200
C   two options for RRR<20 : =10 then linear in log t from 4 to 273
C                            =1
C ***
      ti=dble(tti)
      rrr=dble(rrrr)
      bf=dble(bbf)
      t=log(ti)
      t10=log10(ti)
      rrrl=log(rrr)
      rla=log(1000.)
      rlb=log(500.)
      rlc=log(200.)
      rld=log(100.)
      rle=log(50.)
      rlf=log(20.)
C ** check for rrr<20
      if (rrr.lt.15.)go to 151
      if (rrr.lt.20.)rrrl=rlf
C ** check flat region where r=r273/rrr
      t_flat=4.723436d1-7.201784d00*rrrl+3.075152d-1*rrrl**2
      if(ti.lt.t_flat)then
        rho=r273/rrr
        go to 99
      endif
C ** check to see if t .ge. 273  then rho independent of rrr
151   if(ti.ge.373.)then
         rho=r373+rupp*(ti-373.)
         go to 99
      endif
      if(ti.gt.273.)then
         rho=r273+rup*(ti-273.)
         go to 99
      endif
C  check to see if rrr is a reference curve
C ** if rrr.lt.15 reset to 20, if rrr.gt.1000 reset to 1000
C ** if rrr.lt.15 or lt. 5 use rrr=10 or rrr=1 
      if (rrr.lt.5)then
        rho=exp(r1(t))
        go to 99
      endif
      if(rrr.lt.15.)then
         rho=exp(r10(t))
         go to 99
      endif
      if(rrr.le.20.)then
        rho=exp(r20(t))
        go to 99
      endif
      if(rrr.eq.  50.)then
        rho=exp(r50(t))
        go to 99
      endif
      if(rrr.eq. 100.)then
        rho=exp(r100(t))
        go to 99
      endif
      if(rrr.eq.  200.)then
        rho=exp(r200(t))
        go to 99
      endif
      if(rrr.eq. 500.)then
        rho=exp(r500(t))
        go to 99
      endif
      if(rrr.ge. 1000.)then
        rho=exp(r1000(t))
        go to 99
      endif
C **  find input data for 20<rrr<50
      if(rrr.lt.50.)then
        nd=3
        nc=3
        x(1)=rlf
        x(2)=rle
        x(3)=rld
        y(1)=r20(t)
        y(2)=r50(t)
        y(3)=r100(t)
C ** check for merging
        if(ti.ge.200.)then
          rho=exp(y(2))
          go to 99
         endif
         if(ti.ge.130.)then
           drrr=(rle-rrrl)/(rle-rlf)
           drt=y(1)-y(2)
           rho=exp(y(2)+drrr*drt)
           go to 99
         endif
        call polfit(x,y,nd,nc,ss,b)
        rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2)
        go to 99
      endif
C **  find input data for 50<rrr<100
      if(rrr.lt.100.)then
      fudge=1.-(130.-ti)/50.
        nd=4
        nc=3
        x(1)=rlf
        x(2)=rle
        x(3)=rld
        x(4)=rlc
        y(1)=r20(t)
        y(2)=r50(t)
        y(3)=r100(t)
        y(4)=r200(t)
C ** check for merging
        if(ti.ge.200.)then
          rho=exp(y(2))
          go to 99
         endif
         if(ti.gt.130.)then
        call polfit(x,y,3,2,ss,b)
          rho=exp(b(1)+b(2)*rrrl)
          go to 99
         endif
         if(ti.ge.80.)then
           drrr=(rld-rrrl)/(rld-rle)
           drt=y(2)-y(3)
           rho=exp(y(3)+drrr*drt)
           go to 99
         endif
        call polfit(x,y,nd,nc,ss,b)
        rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2+b(4)*rrrl**3)
        go to 99
      endif
C **  find input data for 100<rrr<200
      if(rrr.lt.200.)then
        nd=4
        nc=3
        x(1)=rle
        x(2)=rld
        x(3)=rlc
        x(4)=rlb
        y(1)=r50(t)
        y(2)=r100(t)
        y(3)=r200(t)
        y(4)=r500(t)
C ** check for merging
        if(ti.ge.130.)then
          rho=exp(y(2))
          go to 99
         endif
         if(ti.ge.80.)then
           drrr=(rlc-rrrl)/(rlc-rld)
           drt=y(2)-y(3)
           rho=exp(y(3)+drrr*drt)
           go to 99
         endif
         if(ti.ge.60.)then
        call polfit(x,y,3,3,ss,b)
          rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2)
          go to 99
         endif
        call polfit(x,y,nd,nc,ss,b)
        rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2+b(4)*rrrl**3)
        go to 99
      endif
C **  find input data for 200<rrr<500
      if(rrr.lt.500.)then
        nd=4
        nc=3
        x(1)=rld
        x(2)=rlc
        x(3)=rlb
        x(4)=rla
        y(1)=r100(t)
        y(2)=r200(t)
        y(3)=r500(t)
        y(4)=r1000(t)
        if(ti.ge.80.)then
          rho=exp(y(2))
          go to 99
         endif
         if(ti.ge.60.)then
           drrr=(rlb-rrrl)/(rlb-rlc)
           drt=y(2)-y(3)
           rho=exp(y(3)+drrr*drt)
           go to 99
         endif
         if(ti.ge.50.)then
        call polfit(x,y,3,3,ss,b)
          rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2)
          go to 99
         endif
        call polfit(x,y,nd,nc,ss,b)
        rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2+b(4)*rrrl**3)
        go to 99
      endif
C **  find input data for 500<rrr<1000
        nd=3
        nc=3
        x(1)=rlc
        x(2)=rlb
        x(3)=rla
        y(1)=r200(t)
        y(2)=r500(t)
        y(3)=r1000(t)
C ** check for merging
        if(ti.ge.60.)then
          rho=exp(y(1))
          go to 99
         endif
         if(ti.ge.50.)then
           drrr=(rla-rrrl)/(rla-rlb)
           drt=y(2)-y(3)
           rho=exp(y(3)+drrr*drt)
           go to 99
         endif
        call polfit(x,y,nd,nc,ss,b)
        rho=exp(b(1)+b(2)*rrrl+b(3)*rrrl**2)
99     rhoc=rho*1.d-6
C  ** magnetoresistance if BF*r273/r(t) .ge. 10
      rtest=r273/rho
      bbt=bf*rtest
      if(bbt.gt.10.)then
        s=log10(rtest*bf)
        dr=-2.662d00+3.168d-1*s+6.229d-1*s**2-1.839d-1*s**3
     1 +1.827d-2*s**4
        drr=10.d00**dr
        rhm=rho*(drr+1.)
        rhoc=rhm*1.d-6
      endif
      rcu=sngl(rhoc)
      return
      end
      function r1000(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2,a3,a4,a5
     1/-8.5844d00,8.5187d-1,7.45395d-1,-2.643894d1,1.00904d1/
      data a6,a7,a8
     1/-1.0944d0,-1.689324d1,3.55861d0/
      at=exp(t)
      if(at.le.20.)then
        r1000=a1+a2*t
        return
      endif
      if(at.lt.30.)then
        r1000=a4+a5*t+a6*t**2
        return
      endif
      if(at.lt.60.)then
        r1000=a7+a8*t
        return
      endif
      if(at.lt.273.)then
        r1000=t273(t)
        return
      endif
      if(at.lt.373.)then
      r1000=log(r273+dr*(at-273.))
      return
      endif
      r1000=log(r373+drr*(at-373.))
      return
      end
      function r500(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2,a3,a4,a5
     1/-1.08815d1,1.8026d0,3.04869d0,-1.526306d1,3.1474d0/
      at=exp(t)
      if(at.lt.17.)then
        r500=log(r273/500.)
        return
      endif
      if(at.le.25.)then
        r500=a1+a2*t
        return
      endif
      if(at.lt.60.)then
        r500=a4+a5*t
        return
      endif
      if(at.lt.273.)then
        r500=t273(t)
        return
      endif
      if(at.lt.373.)then
      r500=log(r273+dr*(at-273.))
      return
      endif
      r500=log(r373+drr*(at-373.))
      return
      end
      function r200(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2,a3,a4,a5
     1/-9.6812d0,1.6508d0,2.21586d0,-1.356154d1,2.7606d0/
      at=exp(t)
      if(at.le.22.)then
        r200=-7.9317d0+1.0848d0*t
        return
      endif
      if(at.lt.33.)then
        r200=a1+a2*t
        return
      endif
      if(at.lt.60.)then
        r200=a4+a5*t
        return
      endif
      if(at.lt.273.)then
        r200=t273(t)
        return
      endif
      if(at.lt.373.)then
      r200=log(r273+dr*(at-273.))
      return
      endif
      r200=log(r373+drr*(at-373.))
      return
      end
      function r100(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2,a3,a4,a5
     1/2.332661d1,-1.712748d1,2.690432d0,-1.208846d1,2.39451d0/
      at=exp(t)
      if(at.le.27.)then
        r100=-6.8338d0+8.909d-1*t
        return
      endif
      if(at.lt.40.)then
        r100=a1+a2*t+a3*t**2
        return
      endif
      if(at.lt.60.)then
        r100=a4+a5*t
        return
      endif
      if(at.lt.273.)then
        r100=t273(t)
        return
      endif
      if(at.lt.373.)then
      r100=log(r273+dr*(at-273.))
      return
      endif
      r100=log(r373+drr*(at-373.))
      return
      end
      function r50(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2,a3,a4,a5
     1/-5.61639d0,6.6564d-1,-7.3304d0,1.1696d0,-1.13669d1/
      data a6,a7,a8,a9
     1/2.2638d0,-6.2344d0,9.67535d-1,3.94729d-2/
      at=exp(t)
      if(at.lt.25.)then
        r50=log(r273/50.)
        return
      endif
      if(at.lt.30.)then
        r50=a1+a2*t
        return
      endif
      if(at.le.40.)then
        r50=a3+a4*t
        return
      endif
      if(at.lt.100.)then
        r50=a5+a6*t
        return
      endif
      if(at.lt.273.)then
        r50=a7+a8*t+a9*t**2
        return
      endif
      if(at.lt.373.)then
      r50=log(r273+dr*(at-273.))
      return
      endif
      r50=log(r373+drr*(at-373.))
      return
      end
      function r20(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2,a3,a4,a5
     1/-5.1921d0,7.912d-1,-4.8473d0,6.8983d-1,-8.000694d0/
      data a6,a7,a8,a9
     1/1.54112d0,-1.36799d0,-9.24664d-1,2.23252d-1/
      at=exp(t)
      if(at.lt.41.)then
        r20=a1+a2*t
        return
      endif
      if(at.lt.100.)then
        r20=a5+a6*t
        return
      endif
      if(at.lt.264.)then
        r20=a7+a8*t+a9*t**2
        return
      endif
      if(at.lt.273.)then
        r20=-3.3743+0.68d0*t
        return
      endif
      if(at.lt.373.)then
      r20=log(r273+dr*(at-273.))
      return
      endif
      r20=log(r373+drr*(at-373.))
      return
      end
      function r10(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      data a1,a2
     1/-2.654d0,5.516d-1/
      at=exp(t)
      if(at.lt.273.)then
        r10=a1+a2*t
        return
      endif
      if(at.lt.373.)then
      r10=log(r273+dr*(at-273.))
      return
      endif
      r10=log(r373+drr*(at-373.))
      return
      end
      function r1(t)
      implicit double precision (a-h,o-z)
      data r273,dr,r373,drr/1.553,4.28d-3,1.981,3.8d-3/
      at=exp(t)
      if(at.lt.273.)then
        r1=log(r273)
        return
      endif
      if(at.lt.373.)then
      r1=log(r273+dr*(at-273.))
      return
      endif
      r1=log(r373+drr*(at-373.))
      return
      end
      function t273(t)
      implicit double precision (a-h,o-z)
      data a1,a2,a3
     1/-2.129467d1,6.71069d0,-5.059114d-1/
      t273=a1+a2*t+a3*t**2
      return
      end

	SUBROUTINE POLFIT(X,Y,NPTS,NTERMS,CHISQQ,COEFF)
      	IMPLICIT REAL*8(A-H,O-Z)
c        REAL X,Y,CHISQQ,COEFF
        DIMENSION COEFF(1)
	DIMENSION X(1),Y(1),COEF(25)
	DIMENSION SUMX(49),SUMY(25),ARRAY(25,25)
	LOGICAL TOOBIG
	NMAX = 2*NTERMS - 1
	DO N = 1,NMAX
	  SUMX(N) = 0
	END DO
	DO J = 1,NTERMS
	  SUMY(J) = 0
	  COEF(J) = 0
	END DO
	CHISQ = 0
	DO I = 1,NPTS
	  XI = X(I)
	  YI = Y(I)
      	WEIGHT = 1.0
	  XTERM = WEIGHT
	  DO N = 1,NMAX
	    SUMX(N) = SUMX(N) + XTERM
	    XTERM = XTERM*XI
	  END DO
	  YTERM = WEIGHT*YI
	  DO N = 1,NTERMS
	    SUMY(N) = SUMY(N) + YTERM
	    YTERM = YTERM*XI
	  END DO
	  CHISQ = CHISQ + WEIGHT*YI**2
	END DO
	DO J = 1,NTERMS
	  DO K = 1,NTERMS
	    N = J + K - 1
	    ARRAY(J,K) = SUMX(N)
	  END DO
	END DO
	CALL DETERM(ARRAY,NTERMS,DELTA,TOOBIG)
	IF (TOOBIG) GO TO 500
	IF (DELTA.EQ.0.) THEN
	  DO J = 1,NTERMS
	    COEF(J) = 0
	  END DO
	  GO TO 999
	END IF
	DO L = 1,NTERMS
	  DO J = 1,NTERMS
	    DO K = 1,NTERMS
	      N = J + K - 1
	      ARRAY(J,K) = SUMX(N)
	    END DO
	    ARRAY(J,L) = SUMY(J)
	  END DO
	  CALL DETERM(ARRAY,NTERMS,DET,TOOBIG)
	  IF (TOOBIG) GO TO 500
	  COEF(L) = DET/DELTA
      END DO
	GO TO 800
500	CONTINUE
	DO J = 1,NTERMS
	  COEF(J) = 0
	END DO
	DO J = 1,NTERMS
	  DO K = 1,NTERMS
	    N = J + K - 1
	    ARRAY(J,K) = SUMX(N)
	  END DO
	END DO
	CALL MATINV(ARRAY,NTERMS,DELTA)
	IF (DELTA.EQ.0.) THEN
	  GO TO 999
	END IF
	DO I = 1,NTERMS
	  DO J = 1,NTERMS
	    COEF(I) = COEF(I) + ARRAY(I,J)*SUMY(J)
	  END DO
	END DO
800	CONTINUE
	DO J = 1,NTERMS
	  CHISQ = CHISQ - 2.D0*COEF(J)*SUMY(J)
	  DO K = 1,NTERMS
	    N = J + K - 1
	    CHISQ = CHISQ + COEF(J)*COEF(K)*SUMX(N)
	  END DO
	END DO
999	CONTINUE
        DO 1 I=1,NTERMS
1       COEFF(I)=COEF(I)
        CHISQQ=CHISQ
        RETURN
      	END



	SUBROUTINE DETERM (ARRAY,NORDER,DET,TOOBIG)
      	IMPLICIT REAL*8(A-H,O-Z)
	REAL*8 ARRAY(25,25),SAVE,DET
	LOGICAL TOOBIG
	TOOBIG = .FALSE.
	DET = 1
	DO K = 1,NORDER
	  IF (ARRAY(K,K).EQ.0) THEN
	    DO J = K,NORDER
	      IF (ARRAY(K,J).NE.0) GO TO 10
	    END DO
	    DET = 0
	    GO TO 999
10	    DO I = K,NORDER
	      SAVE = ARRAY(I,J)
	      ARRAY(I,J) = ARRAY(I,K)
	      ARRAY(I,K) = SAVE
	    END DO
	    DET = DET * -1**(J-K)
	  END IF
	  IF ((DLOG10(DABS(DET))+DLOG10(DABS(ARRAY(K,K)))).LT.38) THEN
	    DET = DET * ARRAY(K,K)
	  ELSE
	    TOOBIG = .TRUE.
	    GO TO 999
	  END IF
	  K1 = K+1
	  DO I = K1,NORDER
	    DO J = K1,NORDER
	      ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)/ARRAY(K,K)*ARRAY(K,J)
	    END DO
	  END DO
	END DO
999	RETURN
	END


	SUBROUTINE MATINV(ARRAY,NORDER,DET)
      	IMPLICIT REAL*8(A-H,O-Z)
	REAL*8 ARRAY(25,25),DET,SAVE,AMAX
	DIMENSION IK(25),JK(25)
	DET = 1
	DO K = 1,NORDER
	  AMAX = 0
	  DO I = K,NORDER
	    DO J = K,NORDER
	      IF (DABS(AMAX).LT.DABS(ARRAY(I,J))) THEN
		AMAX = ARRAY(I,J)
		IK(K) = I
		JK(K) = J
	      END IF
	    END DO
	  END DO
	  IF (AMAX.EQ.0) THEN
	    DET = 0
	    GO TO 999
	  END IF
	  I = IK(K)
	  IF (I.GT.K) THEN
	    DO J = 1,NORDER
	      SAVE = ARRAY(K,J)
	      ARRAY(K,J) = ARRAY(I,J)
	      ARRAY(I,J) = -SAVE
	    END DO
	  END IF
	  J = JK(K)
	  IF (J.GT.K) THEN
	    DO I = 1,NORDER
	      SAVE = ARRAY(I,K)
	      ARRAY(I,K) = ARRAY(I,J)
	      ARRAY(I,J) = -SAVE
	    END DO
	  END IF
	  DO I = 1,NORDER
	    IF (I.NE.K) ARRAY(I,K) = -ARRAY(I,K)/AMAX
	  END DO
	  DO I = 1,NORDER
	    DO J = 1,NORDER
	      IF (I.NE.K.AND.J.NE.K) THEN
		ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K)*ARRAY(K,J)
	      END IF
	    END DO
	  END DO
	  DO J = 1,NORDER
	    IF (J.NE.K) ARRAY(K,J) = ARRAY(K,J)/AMAX
	  END DO
	  ARRAY(K,K) = 1./AMAX
	END DO
	DO L = 1,NORDER
	  K = NORDER - L + 1
	  J = IK(K)
	  IF (J.GT.K) THEN
	    DO I = 1,NORDER
	      SAVE = ARRAY(I,K)
	      ARRAY(I,K) = -ARRAY(I,J)
	      ARRAY(I,J) = SAVE
	    END DO
	  END IF
	  I = JK(K)
	  IF (I.GT.K) THEN
	    DO J = 1,NORDER
	      SAVE = ARRAY(K,J)
	      ARRAY(K,J) = -ARRAY(I,J)
	      ARRAY(I,J) = SAVE
	    END DO
	  END IF
	END DO
999	RETURN
	END