rademit.f

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

	program radiation 
	implicit none
	real*8 aa(10),a,dd,d,pi,w,ww,G,Gint,result
	real*8 bohr,kco,ksi,num,theta,c,bbmin,bmin,fi
	complex*16 eps(10,2),es,krsin,den,Ez0
	real*8 inten0,current,EL,MASS,H,alpha,csi 
	real*8 prefact,probrad,mprobrad,Vb,EPSI0
	real*8 apert,current0,dist,w1,w2,step
	real*8 mreduced
	integer nn,nsurf,m,nmax
	character*1 id,dico
	character*2 sample,tip
	common/sph/a(10),d
	common/omega/ww
	common/surf/nsurf
	common/points/nn(10)
	common/mode/m
	common/diel/eps
	common/reflec/Ez0
	common/hyperb/bmin,fi
	common/bias/Vb
	common/value/result
	common/int/nmax
	data pi/3.14159265358979d0/
	data bohr/0.529177d0/
	data c/137.d0/
	EPSI0=8.854d-12
	EL=1.60219D-19
        MASS=9.10953D-31
        H=1.05459D-34
	csi=2.99792458d8
	alpha=1.d0/c
********************************************************************
**************** Initializing points for integration ***************
********************************************************************
	nmax=64
	call inital(nmax)
********************************************************************
*** Loading the experimental files for the dielectric functions ****
*********** and the Fermi energies of tip and sample ***************
********************************************************************
	call dataexp
********************************************************************
******************** Parameters of the problem *********************
********************************************************************
	nsurf=2
	m=0
	nn(1)=200
	nn(2)=200
	aa(1)=60.d0
	a(1)=aa(1)*10.d0/bohr
	aa(2)=60.d0
	a(2)=aa(2)*10.d0/bohr
	open(unit=3,file='parameters.dat')
	read(3,*)tip
	read(3,*)sample
	call fermien(sample,tip)
	read(3,*) bbmin
	bmin=bbmin/bohr
	read(3,*)apert
	fi=apert*pi/180.d0
	read(3,*)Vb
	read(3,*)id
	dico='d'
	if(id.eq.dico) then
	read(3,*)dist
	dd=dist-bbmin
	d=dd/bohr
	ww=0.d0
	inten0=current(0.d0,0.d0)
	else
	read(3,*)current0
	inten0=current0*1.d-9
	call distance(current0,dist)
	write(6,*)'Current in Amp', inten0,'Distance in Angs',dist
	dd=dist-bbmin
	d=dd/bohr
	endif
	read(3,*)w1
	read(3,*)w2
	read(3,*)step
	close(3)
********************************************************************
************************ Loop of energies **************************
********************************************************************
	open(unit=2,file='output.dat')
	write(2,*)'Current (Amp) and distance (Angs)',inten0,dist
	do w=w1,w2,step
	ww=w/27.2d0
	call dielectric(ww)
********************************************************************
******** Here it is the information about the reflected wave *******
********************************************************************
	theta=45.d0/180.d0*pi
	es=eps(1,2)
	kco=ww/c*dcos(theta)
	ksi=ww/c*dsin(theta)
	krsin=cdsqrt((ww/c)**2*es-ksi**2)
	num=2.d0*es*kco*dsin(theta)
	den=es*kco+krsin
	Ez0=num/den
********************************************************************
******************* Matrix inversion process ***********************
********************************************************************
	call inversion
	call coefficients
	call ginter
********************************************************************
*************** Integration of the local field factor  *************
********************************************************************
	call sumatorio

	Gint=G(dist/bohr)
	
********************************************************************
**************** Calculation of the Radiated Power *****************
********************************************************************
	prefact=((w*EL/H)**2)*EL/((2.d0*pi)**2)/((csi)**3)/H
	probrad=prefact/4.d0/pi/EPSI0*result
	mprobrad=probrad/1.d6
	mreduced=mprobrad/inten0/1.d9
	write(6,*) sngl(w),sngl(Gint),sngl(mprobrad),
     c         sngl(mreduced)
	write(2,*) sngl(w),sngl(Gint),sngl(mprobrad),
     c         sngl(mreduced)
	enddo
	write(2,*)'Energy(eV),Field enhanc, Rad power, Norm rad pow.'
	write(2,*)'Intensity',inten0,' Amps'
	close(2)
	
	end	

********************************************************************
******** Here the field is evaluated *****
********************************************************************
	real*8 function G(z)
	implicit none
	real*8 z,ro,zo,rroo,zzoo,bohr
	complex*16 Ez0,Eziz,pintz,Ezir,pintr
	integer l,i,nn,ntot,nsurf,npar,p
	common/reflec/Ez0
	common/points/nn(10)
	common/surf/nsurf
	common/param/ro,zo
	data bohr/0.529177d0/
********************************************************************
	rroo=0.d0
	ro=rroo*10.d0/bohr
	zzoo=z
	zo=z
********************************************************************
	ntot=0
	do l=1,nsurf
	ntot=ntot+nn(l)
	enddo
	Eziz=(0.d0,0.d0)
	l=1
	npar=nn(l)
	do i=1,ntot
	if (i.eq.npar+1) then
	l=l+1
	npar=0.d0
	do p=1,l
	npar=npar+nn(p)
	enddo
	endif
	call parcialz(l,i,pintz)
	Eziz=Eziz+pintz
	enddo
	Ezir=0.d0
	l=1
	npar=nn(l)
	do i=1,ntot
	if (i.eq.npar+1) then
	l=l+1
	npar=0.d0
	do p=1,l
	npar=npar+nn(p)
	enddo
	endif
	call parcialr(l,i,pintr)
	Ezir=Ezir+pintr
	enddo
	G=dsqrt((cdabs(Ez0+Eziz))**2+(cdabs(Ezir))**2)
	return
	end

******** Here it is calculated the induced field at the point ro,zo 
******** given before ****************************************

	subroutine parcialz(l,i,pintz)
        implicit none
        integer teta,i,nn,nsurf,l,npar,k,nmax,ip,l2
        real*8 r,rprim,zprim,fact1,pi,e,err
        real*8 ro,zo,a,d,fi0,fi1,fintrez,fintimz,fre,gim,tere
        complex*16 coefval,pintz,u,fog
        external fintrez
        external fintimz
	common/param/ro,zo
        common/intco/coefval(0:1000)
        common/sph/a(10),d
        common/points/nn(10)
        common/surf/nsurf
        common/parc/tere
        common/parc2/l2
        data pi/3.14159265358979d0/
        u=(0.d0,1.d0)
	npar=0
        fi0=0.d0
        fi1=2.d0*pi
        err=1.d-6
        nmax=64
        ip=0
        do k=1,l
        npar=npar+nn(k)
        enddo
        teta=i-(npar-nn(l))
        tere=-pi/2.d0+(teta-0.5d0)/nn(l)*pi
        l2=l
        fact1=dsqrt(rprim(l,tere)**2+zprim(l,tere)**2)
        call qfint(fi0,fi1,fintrez,fre,err,e,nmax,ip)
        call qfint(fi0,fi1,fintimz,gim,err,e,nmax,ip)
	fog=dcmplx(fre,gim)
	pintz=r(l,tere)*fact1*fog*coefval(i)*pi/nn(l)
        return
        end
	
	subroutine parcialr(l,i,pintr)
        implicit none
        integer teta,i,nn,nsurf,l,npar,k,nmax,ip,l2
        real*8 r,rprim,zprim,fact1,pi,e,err
        real*8 ro,zo,a,d,fi0,fi1,fintrer,fintimr,fre,gim,tere
        complex*16 coefval,pintr,u,fog
        external fintrer
        external fintimr
	common/param/ro,zo
        common/intco/coefval(0:1000)
        common/sph/a(10),d
        common/points/nn(10)
        common/surf/nsurf
        common/parc/tere
        common/parc2/l2
        data pi/3.14159265358979d0/
        u=(0.d0,1.d0)
	npar=0
        fi0=0.d0
        fi1=2.d0*pi
        err=1.d-6
        nmax=64
        ip=0
        do k=1,l
        npar=npar+nn(k)
        enddo
        teta=i-(npar-nn(l))
        tere=-pi/2.d0+(teta-0.5d0)/nn(l)*pi
        l2=l
        fact1=dsqrt(rprim(l,tere)**2+zprim(l,tere)**2)
        call qfint(fi0,fi1,fintrer,fre,err,e,nmax,ip)
        call qfint(fi0,fi1,fintimr,gim,err,e,nmax,ip)
	fog=dcmplx(fre,gim)
	pintr=r(l,tere)*fact1*fog*coefval(i)*pi/nn(l)
        return
        end

	real*8 function fintrez(fi)
        implicit none
        real*8 fi,fio,r,z,ro,zo,fact1,fact2,teta
        integer l,m
        common/mode/m
        common/parc/teta
        common/parc2/l
	common/param/ro,zo
        fio=0.d0
	fact1=r(l,teta)**2+ro**2-2.d0*r(l,teta)*ro*dcos(fi-fio)
        fact2=(z(l,teta)-zo)**2
        fintrez=-((zo-z(l,teta))/(dsqrt(fact1+fact2))**3)*dcos(m*fi)
        return
        end

	real*8 function fintimz(fi)
        implicit none
        real*8 fi,fio,r,z,ro,zo,fact1,fact2,teta
        integer l,m
        common/mode/m
        common/parc/teta
        common/parc2/l
	common/param/ro,zo
        fio=0.d0
	fact1=r(l,teta)**2+ro**2-2.d0*r(l,teta)*ro*dcos(fi-fio)
        fact2=(z(l,teta)-zo)**2
        fintimz=-((zo-z(l,teta))/(dsqrt(fact1+fact2))**3)*dsin(m*fi)
        return
        end
	
	real*8 function fintrer(fi)
        implicit none
        real*8 fi,fio,r,z,ro,zo,fact1,fact2,teta
        integer l,m
        common/mode/m
        common/parc/teta
        common/parc2/l
	common/param/ro,zo
        fio=0.d0
	fact1=r(l,teta)**2+ro**2-2.d0*r(l,teta)*ro*dcos(fi-fio)
        fact2=(z(l,teta)-zo)**2
        fintrer=-(ro-r(l,teta)*dcos(fi-fio))
     c        /((dsqrt(fact1+fact2))**3)*dcos(m*fi)
        return
        end

	real*8 function fintimr(fi)
        implicit none
        real*8 fi,fio,r,z,ro,zo,fact1,fact2,teta
        integer l,m
        common/mode/m
        common/parc/teta
        common/parc2/l
	common/param/ro,zo
        fio=0.d0
	fact1=r(l,teta)**2+ro**2-2.d0*r(l,teta)*ro*dcos(fi-fio)
        fact2=(z(l,teta)-zo)**2
        fintimr=-(ro-r(l,teta)*dcos(fi-fio))
     c        /((dsqrt(fact1+fact2))**3)*dsin(m*fi)
        return
        end


******** Here the matrix inversion process starts ********
	subroutine inversion
	implicit none
	real*8 pi
	complex*16 Gmatrix,glin,a(1000,1000),c(1000,1),lambda
	integer teta,tetap,i,j,k,l,m,m0,mp,n,ntot,npar,npari,np,nn,nsurf
	common/matrix/glin(0:1000,0:1000)
	common/points/nn(10)
	common/surf/nsurf
	data pi/3.14159265358979d0/
	ntot=0
	do l=1, nsurf
	ntot=ntot+nn(l)
	enddo
	l=1
	npar=nn(1)
	do j=1,ntot
	if(j.eq.npar+1) then
	l=l+1
	npar=0.d0
	do m=1,l
	npar=npar+nn(m)
	enddo
	endif
	k=1
	npari=nn(1)
	do i=1,ntot
	if(i.eq.npari+1) then
	k=k+1
	npari=0.d0
	do m=1,k
	npari=npari+nn(m)
	enddo
	endif
	teta=i-(npari-nn(k))
	tetap=j-(npar-nn(l))
	if(i.eq.j) then
	a(i,i)=lambda(l) - Gmatrix(k,l,teta,tetap)
	else
	a(i,j)= - pi/nn(l)*Gmatrix(k,l,teta,tetap)
	endif
	enddo
	c(j,1)=(1.d0,0.d0)
	enddo
	n=ntot
	np=1000
	mp=1
	m0=1
	call gauss(a,c,n,np,m0,mp)
	do i=1,ntot
	do j=1,ntot
	glin(i,j)=a(i,j)
	enddo
	enddo
	return
	end

******** These are the m-components of the external field that 
******** excite the different modes in the system ************

	subroutine coefficients
	implicit none
	integer teta,i,j,l,nn,nsurf,ntot,npar
	complex*16 glin,coefin,val(0:1000)
	complex*16 coef(0:1000),coefval,indepen
	common/matrix/glin(0:1000,0:1000)
	common/points/nn(10)
	common/surf/nsurf
	common/intco/coefval(0:1000)
	ntot=0
	do l=1,nsurf
	ntot=ntot+nn(l)
	enddo
	l=1
	npar=nn(l)
	do j=1,ntot
	if(j.eq.npar+1) then
	l=l+1
	npar=0.d0
	do i=1,l
	npar=npar+nn(i)
	enddo
	endif
	teta=j-(npar-nn(l))
	coef(j)=indepen(l,teta)
	enddo
	do i=1,ntot
	val(i)=(0.d0,0.d0)
	do j=1,ntot
	coefin=glin(i,j)*coef(j)
	val(i)=val(i)+coefin
	enddo
	coefval(i)=val(i)	
	enddo
	return
	end

	complex*16 function indepen(l,teta)
	implicit none
	real*8 fact1
	real*8 zprim,rprim,teta2,pi
	integer l,teta,l2,nn
	complex*16 fog,Ez0
	common/tirc/l2
	common/tirc2/teta2
	common/points/nn(10)
	common/reflec/Ez0
	data pi/3.14159265358979d0/
	l2=l
	teta2=-pi/2.d0+(teta-0.5d0)/nn(l)*pi
	fact1=dsqrt(rprim(l,teta2)**2+zprim(l,teta2)**2)
	fog=Ez0*(-rprim(l,teta2))
	indepen=fog/fact1
	return
	end


	complex*16 function lambda(l)
	implicit none
	real*8 pi
	integer l
	complex*16 eps
	common/diel/eps(10,2)
	data pi/3.1415926535897932d0/
	lambda=2.d0*pi*(eps(l,2)+eps(l,1))/(eps(l,2)-eps(l,1))
	return
	end