PROGRAM BCS C====================================================================== C THIS PROGRAM CALCULATES THE CONDUCTIVITY (AS A FUNCTION C OF FREQUENCY) FOR A SUPERCONDUCTOR ACCORDING TO THE RAINER- C BRANDT EXTENSION OF THE MATTIS-BARDEEN THEORY. C PHONON LORENTZIANS CAN BE ADDED. C PARAMETERS ARE: TEMPERATURE, CRIT. TEMPERATURE, PLASMAFREQUENCIES C EPSILON INFINITY C NUMBER OF PHONONS, THEIR FREQUENCY, OSCILLATOR-STRENGTH AND C BANDWIDTH C====================================================================== REAL wMIN,wmax,w,EPSINF,gmx,gmn * ,wPLA,wPLA2,wTAU,T,TMN,TSTP,TC,PI,OMEGA COMPLEX EPSIL,S,EPSfc,rp,clad INTEGER ITEMP,NUMPH,I,J,NUM,NTEMP real KT,gapt COMMON /SC/ T,gapt,OMEGA,wTAU,S COMMON /PI/ PI COMMON /ACCURA/ MAX PI=4.*ATAN(1.E0) KT=0 max=8 READ(*,*) T,Tc READ(*,*) EPSINF,wpla,wtau READ(*,*) wmin,wmax,NUM wpla2=wpla**2 gmx=alog10(wmax) gmn=alog10(wmin) do 50 i=1,num w=10**(gmn+(I-1)*(gmx-gmn)/num) * w=wmin+(I-1)*(wmax-wmin)/num IF (T.LT.TC) THEN kt=t gapt=1.76*tc*sqrt(cos(0.5*pi*(T/Tc)**2)) OMEGA=w CALL S1 epsfc=S*(0.,1.)*wPLA2/(wTAU*OMEGA) else epsfc=-wPLA2/(w*CMPLX(w,wTAU)) endif EPSIL=epsinf+EPSfc write(*,*) w,real(epsil),aimag(epsil) 50 continue END c====================================================================== SUBROUTINE S1 c PARAMETER FOR INTEGRATION ROUTINE BERS: c xx=omega/2delta(T), xs=omegatau/2delta(T), t=Temperatur/2delta(T) c c PARAMETERS OBTAINED BY MAIN PROGRAM: c T , Tc, omega, omegat=omegatau c REMARK: THIS SOUBROUTINE DEALS ONLY WITH THE REDUCED c TEMPERATURE,2delta(0) IS FIXED BY THE MAIN PROGRAM, THEREFORE c Tc MUST BE FREE, BECAUSE THIS ROUTINE ASSUMES 2delta(0)/kTc=3.52 c c s=sigma/sigma0 IS RETURNED TO THE MAIN PROGRAM c====================================================================== REAL XX,XS real t,gapt,gap,omega,omegat,tt complex s COMMON /SC/ t,gapt,OMEGA,OMEGAT,S COMMON /PI/ PI COMMON /ACCURA/ M gap=2.*gapt*0.695 c tfac=sqrt(cos(0.5*pi*(tred**2))) c tfac=sqrt(1-tred)*(0.9963 + 0.7733 * tred) tt=t/(2.*gapt) xx = omega/gap xs = omegat/gap call BERS(xx, xs, tt, s) RETURN END c====================================================================== c====================================================================== subroutine BERS(x, xs, t, s) c Berechnet komplexe Leitfaehigkeit s =sigma(q=0, omega) fuer c Supraleiter mit beliebiger Stosszeit taus (London-Fall q=0). c Input: x =omega/2Delta, xs =1/taus2Delta, t=Temperatur/2Delta. c c Da die Energieintegrale (Variable: e=E/2Delta) an den Grenzen c divergieren, wird substituiert e=e(u) mit e'(u)=0 an den Grenzen. c Die neue Variable u laeuft von 0 bis sqrt(emax) oder von 0 bis 1 c mit Schrittweite 1/M oder dx. Z.B.: M=10...30 (sehr genau). c====================================================================== complex s, s1, s2, s3, GK REAL D1,DX,U COMMON /ACCURA/ M d1=1./M dx = 1./int(M*Max(1.,sqrt(x))) s1=(0., 0.) s2=(0., 0.) s3=(0., 0.) do 2 u= dx*.5, 1., dx 2 s2 = s2 +GK(.5 +(u/(1.-u))**2, x, xs, t, 2)*u/(1.-u)**3 s = s2*dx*2. if(x .lt. 1) then do 4 u= dx*.5, 1., dx 4 s1 = s1 +GK(.5 +x*u*u*(3.-u-u), x, xs, t, 1)*u*(1.-u) s = s + s1*dx*6.*x else do 6 u= dx*.5, 1., dx 6 s3 = s3 +GK(.5 +(x-1.)*u*u*(3.-u-u), x, xs, t, 3)*u*(1.-u) do 8 u= d1*.5, 1., d1 8 s1 = s1 +GK(x-.5 +u*u*(3.-u-u), x, xs, t, 1)*u*(1.-u) s = s + (s3*dx*(x-1.) +s1*d1)*6. end if s = s *cmplx(0.,xs)*.5/x end c====================================================================== c********************************************************************** complex function GK(e, x, xs, t, k) c Integranden g1, g2, g3 (=gk, k=1,2,3) complex cs,p4,c42 REAL P1,P2,P3,TH,C12,C32 if(k.eq.2) p1=sqrt( ((e+x)**2 -.25)) p2=sqrt( ( e*e -.25)) if(k.eq.3) p3=sqrt( ((e-x)**2 -.25)) if(k.eq.1) p4=cmplx(0., sqrt( (.25 -(e-x)**2))) cs=cmplx(0., xs) th=tanh(e/(t+t+.001)) if(k.eq.1) c42=(.25 +e*(e-x))/(p4*p2 +1E-20) if(k.eq.2) c12=(.25 +e*(e+x))/(p1*p2 +1E-20) if(k.eq.3) c32=(.25 +e*(e-x))/(p3*p2 +1E-20) if(k.eq.1) GK= th* ((1-c42)/(p4+p2+cs) -(1+c42)/( p4-p2+cs)) if(k.eq.2) GK= tanh((e+x)/(t+t+.001))* & ((1+c12)/(p1-p2+cs) -(1-c12)/(-p1-p2+cs)) & + th* ((1-c12)/(p1+p2+cs) -(1+c12)/( p1-p2+cs)) if(k.eq.3) GK= th* ((1-c32)/(p3+p2+cs) -(1+c32)/( p3-p2+cs)) end C======================================================================