materials.f

计算STM诱导金属表面等离激元发光的程序

	subroutine dielectric(w)
	implicit none
	real*8 w
	integer nsurf
	complex*16 eps(10,2),eps1,eps2
	common/diel/eps
	common/surf/nsurf
	eps(1,1)=(1.d0,0.d0)
	eps(1,2)=eps1(w)
	eps(2,1)=eps2(w)
	eps(2,2)=(1.d0,0.d0)
	return
	end

	subroutine fermien(sample,tip)
	implicit none
	character*2 sample,tip
	character*2 mat1,mat2,mat3,mat4,mat5
	character*2 mat6,mat7
	real*8 EFS,EFT,WS,WT
	common/fermi/EFS,EFT,WS,WT

******  FREE ELECTRON MODEL DATA INPUT
******  USED VALUES ARE  BANDWIDTH (EF)  WORKFUNCTION (W)
******      Au             5.3               4.4
******      W              8.0               5.2
******      Ag             5.5               4.7
******      Ir             9.3               5.8
******      Cu             7.0               4.8
******      Pd             x.x               5.2
******      NiAl          11.2               5.1
******   THE INDEX S REFERS TO THE SAMPLE MATERIAL 
******   INDEX T REFERS TO THE TIP MATERIAL.

	mat1='ag'
	mat2='au'
	mat3='cu'
	mat4='w'
	mat5='ir'
	mat6='nial'
	mat7='pd'
        if(sample.eq.mat1) then
	EFS=5.5d0
	WS=4.7d0
	endif
	if(sample.eq.mat2) then
	EFS=5.3d0
	WS=4.4d0
	endif
	if(sample.eq.mat3) then
	EFS=7.d0
	WS=4.8d0
	endif
	if(sample.eq.mat4) then
	EFS=8.d0
	WS=5.2d0
	endif
	if(sample.eq.mat5) then
	EFS=9.3d0
	WS=5.2d0
	endif
	if(sample.eq.mat6) then
	EFS=11.2d0
	WS=5.1d0
	endif
	if(sample.eq.mat7) then
        EFS=7.6d0
	WS=5.2d0
	endif
	if(tip.eq.mat1) then
	EFT=5.5d0
	WT=4.7d0
	endif
	if(tip.eq.mat2) then
	EFT=5.3d0
	WT=4.4d0
	endif
	if(tip.eq.mat3) then
	EFT=7.d0
	WT=4.8d0
	endif
	if(tip.eq.mat4) then
	EFT=8.d0
	WT=5.2d0
	endif
	if(tip.eq.mat5) then
	EFT=9.3d0
	WT=5.2d0
	endif
	if(tip.eq.mat6) then
	EFT=11.2d0
	WT=5.1d0
	endif
	if(tip.eq.mat7) then
	EFT=7.6d0
	WT=5.2d0
	endif

	return
	
	end

c *******************************************************************
c .................................................................
c .................................................................

	subroutine dataexp
c	Experimental values for both real and imaginary parts of e(w) are readed
c	A splin procedure is used to interpolated values.

	implicit real*8(a-h,o-z)
	
	common/resp11/w1(200),r1(200),y1(200),zn1
	common/resp12/f11(200),f12(200)    

	common/resp21/w2(200),r2(200),y2(200),zn2
	common/resp22/f21(200),f22(200)    

c  Reading the experimental data of Ag. 
******** Ag from Peter Johansson's data ********

        open(unit=1,file='epssample.dat',status='old')

	read(1,*) n1
	zn1=dfloat(n1)
	
	do i=1,n1
	   read(1,*) w1(i),r1(i),y1(i)
	   w1(i)=w1(i)/27.2
	end do
        close(unit=1)

	d0=(r1(2)-r1(1))/(w1(2)-w1(1))
	df=(r1(n1)-r1(n1-1))/(w1(n1)-w1(n1-1))
	call spline(w1,r1,n1,d0,df,f11)

	d0=(y1(2)-y1(1))/(w1(2)-w1(1))
	df=(y1(n1)-y1(n1-1))/(w1(n1)-w1(n1-1))
	call spline(w1,y1,n1,d0,df,f12)

c Reading Ir data. 
******** Ir from Peter Johansson's data ********
      	
	open(unit=3,file='epstip.dat',status='old')

       	read(3,*) n2
	zn2=dfloat(n2)
	
	do i=1,n2
 	  read(3,*) w2(i),r2(i),y2(i)
	   w2(i)=w2(i)/27.2
	end do
        close(unit=3)

	d0=(r2(2)-r2(1))/(w2(2)-w2(1))
	df=(r2(n2)-r2(n2-1))/(w2(n2)-w2(n2-1))
	call spline(w2,r2,n2,d0,df,f21)

	d0=(y2(2)-y2(1))/(w2(2)-w2(1))
	df=(y2(n2)-y2(n2-1))/(w2(n2)-w2(n2-1))
	call spline(w2,y2,n2,d0,df,f22)

	return

c *************************************************************
c *************************************************************
        end
c ........................................................................
c ........................................................................
c ........................................................................
	complex*16 function eps1(z)

	implicit real*8(a-h,o-z)
	
	common/resp11/w(200),r(200),y(200),zn
	common/resp12/f1(200),f2(200)    
	complex u1
	data u1/(0.d0,1.d0)/

	n=int(zn) 

	call splint(w,r,f1,n,z,rf)
	call splint(w,y,f2,n,z,yf)

	eps1=rf+yf*u1
	end

c ........................................................................
	complex*16 function eps2(z)

	implicit real*8(a-h,o-z)
	
	common/resp21/w(200),r(200),y(200),zn
	common/resp22/f1(200),f2(200)    
	complex u1
	data u1/(0.d0,1.d0)/

	n=int(zn) 

	call splint(w,r,f1,n,z,rf)
	call splint(w,y,f2,n,z,yf)

	eps2=rf+yf*u1
	end

c........................................................................
c 	Subroutina para interpolacion por un splin cubico

	subroutine spline(x,y,n,yp1,ypn,y2)
 	implicit real*8 (a-h,o-z)
	parameter (nmax=201)
	dimension x(n),y(n),y2(n),u(nmax)
	if (yp1.gt..99d30) then
		y2(1)=0.d0
		u(1)=0.d0
	else
		y2(1)=-.5d0
		u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
	endif
	do 11 i=2,n-1
		sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
		p=sig*y2(i-1)+2.d0
		y2(i)=(sig-1.d0)/p
		u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
     *		     /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11	continue
	if(ypn.gt..99d30) then
		qn=0.d0
		un=0.d0
	else
	    qn=0.5d0
	    un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
	endif
	y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)

	do 12 k=n-1,1,-1
		y2(k)=y2(k)*y2(k+1)+u(k)
12	continue
	return
	end


		
	subroutine splint(xa,ya,y2a,n,x,y)
	implicit real*8(a-h,o-z)
	dimension xa(n), ya(n), y2a(n)
	klo=1
	khi=n
1	if (khi-klo.gt.1) then
		k=(khi+klo)/2
		if(xa(k).gt.x) then
			khi=k
		else
			klo=k
		endif
		goto 1
	endif
	h=xa(khi)-xa(klo)
	if(h.eq.0) pause 'bad xa input'
	a=(xa(khi)-x)/h
	b=(x-xa(klo))/h
	y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))
     *  *(h**2)/6
	return
     	end