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