c ###################################################################### c # c # NB3SN PROPERTIES PACKAGE - ITER SCALING c # --------------------------------------- c # Contains functions for the calculation of the critical properties c # and thermal properties of Nb3Sn c # c ###################################################################### c ###################################################################### block data Nb3SnConstants c ###################################################################### c # c # Niobium-Tin (Nb3Sn) - ITER Jc parameterization c # c # The constants are based on a fit of the ITER CS wire adjusted to c # yield Jc=1000 A/mm**2 at 4.2 K, 12 T and 0 % intrinsic strain c # c # real Tc0Nb3Sn,Bc20Nb3Sn,C0Nb3Sn real pNb3Sn,qNb3Sn,nuNb3Sn real Ca1Nb3Sn,Ca2Nb3Sn,emaxNb3Sn,e0aNb3Sn c * common /Nb3SnData/ Tc0Nb3Sn,Bc20Nb3Sn,C0Nb3Sn, & pNb3Sn,qNb3Sn,nuNb3Sn, & Ca1Nb3Sn,Ca2Nb3Sn,emaxNb3Sn,e0aNb3Sn c * data Tc0Nb3Sn /16.73/ data Bc20Nb3Sn/30.23/ data C0Nb3Sn /65.26e9/ data pNb3Sn / 0.56/ data qNb3Sn / 1.75/ data nuNb3Sn / 1.52/ data Ca1Nb3Sn /45.46/ data Ca2Nb3Sn / 6.52/ data emaxNb3Sn/ 0.0/ data e0aNb3Sn / 0.00244/ c * end c ###################################################################### real function dNb3Sn(T) c ###################################################################### c # c # Density of Nb3Sn (stoichiometric) c # c # Range: 0 <= T <= inf K c # c # References c # ---------- c # Y. Iwasa, Case Studies in Superconducting Magnets: Design and c # Operational Issues, Springer c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # T x absolute temperature K c # dNb3Sn x density Kg/m**3 c # c # c # Author : L.Bottura at Cryosoft c # Version: 1.1 November 2013 c # c ###################################################################### implicit none c * external variables real T c * fit variables c * local variables c * dNb3Sn = 8040.0 c * return end c ###################################################################### real function cNb3Sn(T,B,Epslon, & Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c ###################################################################### c # c # Specific Heat of binary Nb3Sn in J/Kg K as a function of temperature c # field and strain for 2 K <= T <= 400 K. The normal conducting c # specific heat and the superconducting jump at the critical c # temperature are obtained from fit data, and in superconducting state c # a simple T**3 dependence is assumed, as measured by Khlopin. No c # allowance is made for mixed state. c # c # Range: 2 <= T <= 400 K c # c # References c # ---------- c # M.N. Khlopin, The Specific Heat of Nb3Sn in Magnetic Fields up to c # 19 T, Sov. Phys. JETP.63.1.164.1986 c # G.S. Knapp, S.D. Bader, Z. Fisk, Phonon Properties of A-15 c # Superconductors Obtained from Heat-Capacity Measurements, Phys. c # Rev. B, 13(9), 3783, 1976. c # MATPRO, quoting: H. Brechna, "Superconducting Magnet Systems, c # Springer, 1979 (2K < T < 28 K); Bubble Chamber Group Data Handbook - c # Selected Cryogenic Data Notebook, J.E. Jensen, R.B. Stewart, W.A. c # Tuttle, BNL, 1966 (28 K < T < 100 K) c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # T x temperature K c # B x magnetic field T c # Epslon x strain - c # Tc0m x critical temperature (b=0,e=0) K c # Bc20m x upper critical field (t=0,e=0) T c # nu x fitting constant for h(t) - c # Ca1 x fitting constant for s(e) - c # Ca2 x fitting constant for s(e) - c # emax x fitting constant for s(e) - c # e0a x fitting constant for s(e) - c # cNb3Sn x specific heat J/kg K c # c # Other functions called: TcNb3Sn,BcNb3Sn c # c # Author : L.Bottura at CryoSoft c # Version: 4.1 November 2013 c # c ###################################################################### implicit none c * real T,B,Epslon c * real Tc0m,Bc20m real nu,Ca1,Ca2,emax,e0a c * real TcNb3Sn,BcNb3Sn c * fit variables real A1,A3,A5 real AA,CC,DD,a,c,d,na,nc,nd,T0 real C1,n1,C2,n2 real Tmin,Tmax data A1 / 0.111323603 / , A3 / 0.001560141 / , & A5 /-6.54253E-07 / data AA / 1531.184988 / , & CC /-1222.227007 / , DD / -1.269824951 / data a / 13.01639235 / , & c / 4.534367264 / , d / 72.12056165 / data na / 2.38482719 / , & nc / 4.874967323 / , nd / 3.98653929 / data T0 / 17.02416462 / data C1 / 0.12 / , C2 / 0.16 / data n1 / 1 / , n2 / 0.33 / data Tmin/ 2.0/, Tmax / 400.0/ save c * local variables real Cn,Cs,dC,Tc,Bc20,Cp real TT,bb c * TT=T TT=min(TT,Tmax) TT=max(TT,Tmin) c * critical temperature Tc = TcNb3Sn(B,Epslon,Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) if(TT.le.Tc) then c * superconducting state Bc20 = BcNb3Sn(0.0,Epslon,Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) bb = B/Bc20 c * compute normal conducting specific heat at Tc if (Tc .le. T0) then Cn= A1*Tc + A3*Tc**3 + A5*Tc**5 else Cn = AA*Tc**na/(a+Tc)**na + & CC*Tc**nc/(c+Tc)**nc + DD*Tc**nd/(d+Tc)**nd endif c * compute jump in specific heat at Tc dC = (C1*(1-bb**n1) + C2*(1-bb**n2)) * Tc c * compute maximum value of superconducting specific heat Cs = Cn+dC c * actual specific heat, scaled to the temperature Cp = Cs * (TT/Tc)**3 else c * normal state if (TT .le. T0) then Cn= A1*TT + A3*TT**3 + A5*TT**5 else Cn = AA*TT**na/(a+TT)**na + & CC*TT**nc/(c+TT)**nc + DD*TT**nd/(d+TT)**nd endif c * Cp = Cn endif c * final value cNb3Sn = Cp c * return end c ###################################################################### real function kNb3Sn(T) c ###################################################################### c # c # Thermal conductivity of Nb3Sn as a function of temperature T, for c # 4 <= T <= 300 K. NOTE: the conductivity is assumed to be continuous at c # the superconducting to normal state transition, mainly because of c # lack of directly measured data. c # c # Range: 2 <= T <= 300 K c # c # References c # ---------- c # P. Bauer, H. Rajainmaki, E. Salpietro, EFDA Material Dat Compilation c # for Superconductor Simulation, EFDA CSU, Garching, 18 April 2007 c # G.E. Childs, L.J. Ericks, R.L Powell, NBS Monograph 131, 1973 c # G.D. Cody, R.W. Cohen, Thermal Condctivity of Nb3Sn, Rev. Mod. c # Phys., 1964 c # H. Brechna, Superconducting Magnets Systems, 1973 c # MATPRO, quoting H. Brechna, Superconducting Magnets Systems, 1973 c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # T x absolute temperature K c # kNb3Sn x thermal conductivity W/m K c # c # Author : L.Bottura & C.Rosso at Cryosoft c # Version: 4.1 November 2013 c # c ###################################################################### implicit none real T c * fit variables real Tmin,Tmax real k0,k1,k2,km real Tm real alpha,beta,n,m data k0 / 0.005692381 / , k1 / 0.019093063 /, & k2 /-0.001795443 / , km / 2.046939922 / data Tm / 21.70346889 / data alpha / 1.864885217 / , beta / 0.625548295 /, & n / 0.317824884 / , m / 3.446828156 / data Tmin/ 2.0/, Tmax / 300.0/ save c * local variables real TT c * TT=T TT=min(TT,Tmax) TT=max(TT,Tmin) c * kNb3Sn = (k0+k1*(TT/Tm)+k2*(TT/Tm)**2) + & 3.0 * km / (alpha*(TT/Tm)**n + beta/(TT/Tm)**m) c * return end c ###################################################################### real function rNb3Sn(T) c ###################################################################### c # c # Resistivity of Nb3Sn c # c # Range: 20 <= T <= 900 K c # c # References c # ---------- c # c # Z. Fisk, G.W. Webb, Saturation of the High-temperature Normal-State c # Resistivity of Superconductors, Phys. Rev. Lett., 36, 1084, 1976 c # compiled by P. Bauer, H. Rajainmaki, E. Salpietro, EFDA Material c # Data Compilation for Superconductor Simulation, EFDA CSU, Garching, c # 18 April 2007 c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # T x absolute temperature K c # rNb3Sn x resistivity Ohm m c # c # Author : L.Bottura at Cryosoft c # Version: 1.0 March 2013 c # c ###################################################################### implicit none c * external variables real T c * fit variables real AA,a,na,B,C real Tmin,Tmax data AA / 8.506987E-7 / data a / 32.60893849 / data na / 2.680037713 / data B / 3.72872E-08 / , C / 4.49E-10 / data Tmin/ 20.0 /, Tmax / 900.0 / c * local variables real TT c * TT=T TT=min(TT,Tmax) TT=max(TT,Tmin) c * rNb3Sn = AA * TT**na / (a+TT)**na + B + C*TT c * return end c ###################################################################### real function TcNb3Sn(B,Epslon, & Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c ###################################################################### c # c # Critical temperature, in K, for Nb3Sn as a function of applied field c # and strain. The Nb3Sn material is mainly characterized by the c # parameters Tc0m and Bc20m. The additional constants are used for the c # details of the scaling function. c # c # References c # ---------- c # L. Bottura, B. Bordini, Jc(B,T,e) Parameterization for the ITER c # Nb3Sn Production, IEEE Trans. Appl. Sup., 19(2), 1477-1480, 2009 c # c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # B x magnetic field T c # Epslon x strain - c # Tc0m x critical temperature (b=0,e=0) K c # Bc20m x upper critical field (t=0,e=0) T c # nu x fitting constant for h(t) - c # Ca1 x fitting constant for s(e) - c # Ca2 x fitting constant for s(e) - c # emax x fitting constant for s(e) - c # e0a x fitting constant for s(e) - c # TcNb3Sn x critical temperature K c # c # Other functions called: sNb3Sn c # c # Author : L. Bottura at CryoSoft c # Version: 1.0 November 2012 c # c ###################################################################### implicit none c * real B,Epslon c * real Tc0m,Bc20m real nu,Ca1,Ca2,emax,e0a c * real Blim,ee,s,Bc20,bb0,tt,Tc real Blow parameter(Blow=1.0e-3) c * real sNb3Sn c * set the lower limit for the field Blim=amax1(abs(B),Blow) c * compute s(e) ee = Epslon-emax s = sNb3Sn(ee,Ca1,Ca2,e0a) c * compute the critical field at zero temperature Bc20 = Bc20m * s c * compute the reduced field bb0 = Blim/Bc20 c * compute the reduced critical temperature if(bb0.gt.0.0 .and. bb0.lt.1.0) then tt = (1-bb0)**(1.0/nu) elseif(bb0.eq.0.0) then tt = 1.0 elseif(bb0.gt.1.0) then tt = 0.0 endif c * compute Tc Tc = Tc0m*s**0.333333333 * tt c * assign the value TcNb3Sn = Tc c * return end c ###################################################################### real function BcNb3Sn(T,Epslon, & Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c ###################################################################### c # c # Critical field, in T, for Nb3Sn as a function of temperature and c # strain. The Nb3Sn material is mainly characterized by the parameters c # Tc0m and Bc20m. The additional constants are used for the details of c # the scaling function. c # c # References c # ---------- c # L. Bottura, B. Bordini, Jc(B,T,e) Parameterization for the ITER c # Nb3Sn Production, IEEE Trans. Appl. Sup., 19(2), 1477-1480, 2009 c # c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # T x temperature K c # Epslon x strain - c # Tc0m x critical temperature (b=0,e=0) K c # Bc20m x upper critical field (t=0,e=0) T c # nu x fitting constant for h(t) - c # Ca1 x fitting constant for s(e) - c # Ca2 x fitting constant for s(e) - c # emax x fitting constant for s(e) - c # e0a x fitting constant for s(e) - c # BcNb3Sn x critical field T c # c # Other functions called: sNb3Sn c # c # Author : L.Bottura at Cryosoft c # Version: 0.0 26.3.2010 c # c ###################################################################### implicit none c * real T,Epslon c * real Tc0m,Bc20m real nu,Ca1,Ca2,emax,e0a c * real Tlim,ee,s,Tc0,tt,bb,Bc real Tlow parameter(Tlow=0.0) c * real sNb3Sn c * set the lower limit for the field Tlim=amax1(T,Tlow) c * compute s(e) ee = Epslon-emax s = sNb3Sn(ee,Ca1,Ca2,e0a) c * compute the critical temperature at zero background field Tc0 = Tc0m * s**0.33333333 c * compute the reduced temperature tt = Tlim/Tc0 c * compute the reduced critical field if(tt.gt.0.0 .and. tt.lt.1.0) then bb = (1-tt**nu) elseif(tt.le.0.0) then bb = 1.0 elseif(tt.gt.1.0) then bb = 0.0 endif c * compute Bc Bc = Bc20m*s * bb c * assign the value BcNb3Sn = Bc c * return end c ###################################################################### real function JcNb3Sn(T,B,Epslon, & Tc0m,Bc20m,C0,p,q,nu,Ca1,Ca2,emax,e0a) c ###################################################################### c # c # Critical (non-copper) current density, in A/m**2, for Nb3Sn as a c # function of temperature, field and strain. The Nb3Sn material is c # mainly characterized by the parameters Tc0m and Bc20m. The constant c # C0 determines the overall scaling of the Jc curve. The additional c # constants are used for the details of the scaling function. c # c # References c # ---------- c # L. Bottura, B. Bordini, Jc(B,T,e) Parameterization for the ITER c # Nb3Sn Production, IEEE Trans. Appl. Sup., 19(2), 1477-1480, 2009 c # c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # T x absolute temperature K c # B x magnetic field T c # Epslon x strain - c # Tc0m x critical temperature (b=0,e=0) K c # Bc20m x upper critical field (t=0,e=0) T c # C0 x normalization constant A T/m**2 c # p x fitting constant for fp(b) - c # q x fitting constant for fp(b) - c # nu x fitting constant for h(t) - c # Ca1 x fitting constant for s(e) - c # Ca2 x fitting constant for s(e) - c # emax x fitting constant for s(e) - c # e0a x fitting constant for s(e) - c # JcNb3Sn x critical current density A/m**2 c # c # Other functions called: TcNb3Sn,BcNb3Sn,sNb3Sn,hNb3Sn,fpNb3Sn c # c # Author : L.Bottura at Cryosoft c # Version: 0.0 26.3.2010 c # c ###################################################################### implicit none c * real T,B,Epslon c * real Tc0m,Bc20m,C0 real p,q,nu,Ca1,Ca2,emax,e0a c * real Blim,Tlim,ee,s,Tc0,Tc,tt,Bc,bb,h,fp,Jc real Blow parameter(Blow=1.0e-3) real Tlow parameter(Tlow=0.0) c * real sNb3Sn,hNb3Sn,fpNb3Sn,TcNb3Sn,BcNb3Sn c * set the lower limit for the field Blim=amax1(abs(B),Blow) c * set the lower limit for the temperature Tlim=amax1(T,Tlow) c * compute s(e) ee = Epslon-emax s = sNb3Sn(ee,Ca1,Ca2,e0a) c * compute the critical temperature at zero background field Tc0 = TcNb3Sn(0.0,Epslon,Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c * compute the reduced temperature tt = Tlim/Tc0 c * trap case of temperature in excess of the critical temperature if(tt.gt.1.0) then JcNb3Sn = 0.0 return endif c * compute the critical field at the temperature and strain Bc = BcNb3Sn(T,Epslon,Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c * compute the reduced field bb = Blim/Bc c * trap case of field in excess of the critical field if(bb.gt.1.0) then JcNb3Sn = 0.0 return endif c * compute h(t) h = hNb3Sn(tt,nu) c * compute fp(b) fp = fpNb3Sn(bb,p,q) c * compute Jc Jc = C0/Blim * s * h * fp c * assign the value JcNb3Sn = Jc c * return end c ###################################################################### real function TcsNb3Sn(B,Epslon,Jop, & Tc0m,Bc20m,C0,p,q,nu,Ca1,Ca2,emax,e0a) c ###################################################################### c # c # Current sharing temperature, in K, for Nb3Sn as a function of field, c # strain, and operating current density. and strain. The Nb3Sn c # material is mainly characterized by the parameters Tc0m and Bc20m. c # The constant C0 determines the overall scaling of the Jc curve. The c # additional constants are used for the details of the scaling c # function. This function is the inverse of the Jc function JcNb3Sn c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # B x magnetic field T c # Epslon x strain - c # Jop x operating current density A/m**2 c # Tc0m x critical temperature (b=0,e=0) K c # Bc20m x upper critical field (t=0,e=0) T c # C0 x normalization constant A T/m**2 c # p x fitting constant for fp(b) - c # q x fitting constant for fp(b) - c # nu x fitting constant for h(t) - c # Ca1 x fitting constant for s(e) - c # Ca2 x fitting constant for s(e) - c # emax x fitting constant for s(e) - c # e0a x fitting constant for s(e) - c # TcsNb3Sn x current sharing temperature K c # c # Other functions called: TcNb3Sn,BcNb3Sn,JcNb3Sn c # c # Author : L.Bottura at Cryosoft c # Version: 0.0 26.3.2010 c # c ###################################################################### implicit none c * real B,Epslon,Jop c * real Tc0m,Bc20m,C0 real p,q,nu,Ca1,Ca2,emax,e0a c * real Bc20,T,Jc,Tc0,tt,ttlow,ttup real error,tolerance data tolerance/1.0e-5/ logical converged c * real TcNb3Sn,BcNb3Sn,JcNb3Sn c ---------------------------------------------------------------------- c write(6,*) 'B ',B c write(6,*) 'Epslon ',Epslon c write(6,*) 'Jop ',Jop c write(6,*) 'Tc0m ',Tc0m c write(6,*) 'Bc20m ',Bc20m c write(6,*) 'C0 ',C0 c write(6,*) 'p ',p c write(6,*) 'q ',q c write(6,*) 'nu ',nu c write(6,*) 'Ca1 ',Ca1 c write(6,*) 'Ca2 ',Ca2 c write(6,*) 'emax ',emax c write(6,*) 'e0a ',e0a c ---------------------------------------------------------------------- c * compute critical field and zero temperature at applied strain T = 0.0 Bc20 = BcNb3Sn(T,Epslon,Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c ---------------------------------------------------------------------- c write(6,*) 'TcsNb3Sn: Bc20',Bc20 c ---------------------------------------------------------------------- c * check that the field is below the upper critical value if(B.ge.Bc20) then TcsNb3Sn=0.0 return endif c * compute critical current density at applied field, strain and c * zero temperature T = 0.0 Jc = JcNb3Sn(T,B,Epslon,Tc0m,Bc20m,C0,p,q,nu,Ca1,Ca2,emax,e0a) c ---------------------------------------------------------------------- c write(6,*) 'TcsNb3Sn: Jc0',Jc c ---------------------------------------------------------------------- c * check that Jop is below the maximum possible Jc if(Jop.ge.Jc) then TcsNb3Sn=0.0 return endif c * compute critical temperature at applied strain for the maximum c * value of Tcs (Tc) Tc0 = TcNb3Sn(B,Epslon,Tc0m,Bc20m,nu,Ca1,Ca2,emax,e0a) c ---------------------------------------------------------------------- c write(6,*) 'TcsNb3Sn: Tc0',Tc0 c ---------------------------------------------------------------------- c * trap for the jop=0 case to avoid iteration problems if(Jop.le.0.0) then TcsNb3Sn = Tc0 return endif c * upper and lower limits for the iterative root search ttup =Tc0/Tc0m ttlow=0.0 c * find the normalised temperature tcs/tc0 by iteration converged=.false. do while(.not.converged) c * root search by straight bisection (most stable way) tt=0.5*(ttlow+ttup) c ---------------------------------------------------------------------- c write(6,*) 'TcsNb3Sn: tt',tt c ---------------------------------------------------------------------- c * compute critical current T = tt * Tc0m Jc = JcNb3Sn(T,B,Epslon,Tc0m,Bc20m,C0,p,q,nu,Ca1,Ca2,emax,e0a) c * decision process if(Jc.gt.Jop) then ttlow = tt elseif(Jc.le.Jop) then ttup = tt elseif(Jc.eq.Jop) then ttup = tt ttlow = tt endif c * check convergence error = abs(ttup-ttlow) converged = error.le.tolerance enddo c * convert normalised to absolute temperature TcsNb3Sn = tt*Tc0m c ---------------------------------------------------------------------- c write(6,*) 'TcsNb3Sn: Tcs',TcsNb3Sn c ---------------------------------------------------------------------- c * return end c ###################################################################### real function fpNb3Sn(b,p,q) c ###################################################################### c # c # Function fp(b) for Nb3Sn, as from the ITER-2008 parameterization c # c # References c # ---------- c # L. Bottura, B. Bordini, Jc(B,T,e) Parameterization for the ITER c # Nb3Sn Production, IEEE Trans. Appl. Sup., 19(2), 1477-1480, 2009 c # c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # b x reduced field - c # p x fitting constant - c # q x fitting constant - c # fpNb3Sn x function fp(b) - c # c # Other functions called: NONE c # c # Author : L.Bottura at Cryosoft c # Version: 0.0 26.3.2010 c # c ###################################################################### implicit none c * real b real p,q c * real fp c * compute fp(b) if(b.gt.0.0 .and. b.lt.1.0) then fp = b**p * (1-b)**q else fp = 0.0 endif c * assign the value fpNb3Sn = fp c * return end c ###################################################################### real function hNb3Sn(t,nu) c ###################################################################### c # c # Function h(t) for Nb3Sn, as from the ITER-2008 parameterization c # c # References c # ---------- c # L. Bottura, B. Bordini, Jc(B,T,e) Parameterization for the ITER c # Nb3Sn Production, IEEE Trans. Appl. Sup., 19(2), 1477-1480, 2009 c # c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # t x reduced temperature - c # nu x fitting constant - c # hNb3Sn x function h(t) - c # c # Other functions called: NONE c # c # Author : L.Bottura at Cryosoft c # Version: 0.0 26.3.2010 c # c ###################################################################### implicit none c * real t real nu c * real h c * compute h(t) if(t.gt.0.0 .and. t.lt.1.0) then h = (1.0-t**nu) * (1-t*t) elseif(t.le.0.0) then h = 1.0 elseif(t.ge.1.0) then h = 0.0 endif c * assign the value hNb3Sn = h c * return end c ###################################################################### real function sNb3Sn(e,Ca1,Ca2,e0a) c ###################################################################### c # c # Function s(eps) for Nb3Sn, as from the ITER-2008 parameterization c # c # References c # ---------- c # L. Bottura, B. Bordini, Jc(B,T,e) Parameterization for the ITER c # Nb3Sn Production, IEEE Trans. Appl. Sup., 19(2), 1477-1480, 2009 c # c # c # variable I/O meaning units c # -------------------------------------------------------------------- c # e x reduced strain - c # Ca1 x fitting constant - c # Ca2 x fitting constant - c # e0a x fitting constant - c # sNb3Sn x function s(eps) - c # c # Other functions called: NONE c # c # Author : L.Bottura at Cryosoft c # Version: 0.0 26.3.2010 c # c ###################################################################### implicit none c * real e real Ca1,Ca2,e0a c * real esh,s,esh2,e0a2,ed,ed2 c * compute s(eps) esh = ca2 * e0a / sqrt(ca1*ca1-ca2*ca2) ed = e-esh esh2 = esh*esh e0a2 = e0a*e0a ed2 = ed*ed s = 1.0 + (ca1 * (sqrt(esh2+e0a2) - sqrt(ed2+e0a2)) - Ca2*e)/ & (1 - Ca1*e0a) c * limit the value of s to zero if(s.lt.0.0) s = 0.0 c * assign the value sNb3Sn = s c * return end