C	(C) Copr. 06/2005  Alexander Moroz
C	Wave-scattering.com
C	http://www.wave-scattering.com/
C	wavescattering@yahoo.com
C
C	Last update: February 2005
C
C	This is EFFE2P, a computer code to calculate the effective
C       medium parameters within the context of the extended
C       Maxwell - Garnett theory (a quasi-static extention of the 
c       Maxwell Garnett theory) in the case of a binary composite
C       (two different types of spheres embedded in a matrix)
C
C       This code has been employed in the article
C
C	V. Yannopapas and A. Moroz, Negative refractive index metamaterials 
C       from inherently non-magnetic materials for deep infrared to terahertz 
C       frequency ranges, published in
C        J. Phys.: Condens. Matter. <b>17</b>, 3717-3734 (2005) 
C       (http://www.iop.org/EJ/abstract/0953-8984/17/25/002).
C
C       where you can also find an additional teoretical background for
C       the extended Maxwell - Garnett theory employed here. 
C
C  The program default options correspond to the cases considered in the above
C  article. (The generated results are shown in Fig. 4a of the article.)
C
C  In order to modify the default option, read through the comments
C  in the program and description of the parameters and variables.
C 
C  OUTPUT:
C  Returns file effmed.dat containing the effective medium parameters 
C  versus wavelength.
C =====
C
C  TESTS:
C
C  1) THE RESULTS OBTAINED BY DEFAULT SETUP OF PARAMETERS
C     SHOULD COINCIDE WITH THOSE DISPLAYED IN FIG. 4a of the article.
C     In particular, the first result lines of effmed.dat with 
C     the program default options should read as follows:
C 
C   13209.000      -50.1126     80.2440        1.3153   0.00080        5.4064      9.7573
C   13206.646      -51.6593     69.7463        1.3162   0.00080        4.8056      9.5469
C   13204.292      -51.5292     60.5496        1.3171   0.00081        4.2896      9.2908
C   13201.939      -50.3046     52.6734        1.3180   0.00081        3.8505      9.0095
C   13199.587      -48.4108     46.0136        1.3189   0.00082        3.4787      8.7172
C   13197.236      -46.1410     40.4189        1.3198   0.00082        3.1645      8.4229
C   13194.886      -43.6896     35.7316        1.3208   0.00083        2.8993      8.1326
C   13192.536      -41.1812     31.8070        1.3217   0.00083        2.6756      7.8495
C   13190.188      -38.6937     28.5202        1.3226   0.00084        2.4876      7.5756
C   13187.840      -36.2743     25.7673        1.3236   0.00084        2.3301      7.3118
C
C========================================================================
C
C   	REQUIREMENTS & NOT DISTRIBUTED DEPENDENT ROUTINES:
C
C       You should supply your own routine
C
C	GNZBESS(RX(I),LMAX,jl,drjl,nl,drnl)
C
C	which generates the Bessel functions JL(0:LMAX),NL(0:LMAX)
C	of the first and second order and their derivatives
C	DJL(0:LMAX),DNL(0:LMAX) for a complex argument RX(I) in the
C       subroutine SPHEREF. (SPHEREF routine calculates the single 
C       sphere scattering properties (including coated spheres) and
C       returns the elements of the diagonal T matrix.)
C	
C	My program was tested together with the AMOS package of 
C	complex Bessel functions routines which 
C	you can, together with their dependencies,
C	freely download at Netlib (http://www.netlib.org/amos).
C	Before using the AMOS package, machine constants
C	in the real function D1MACH(I) and in the integer function I1MACH(I)
C	have to be adjusted to your computer. When using Microsoft Windows,
C	you can activate the same machine constants as listed therein under
C	the option of Silicon Graphics IRIS. 
C
C =====================================================================
C
C	This code is made available on the following conditions:
C
C	- You can use and modify it, as long as acknowledgment is given
C	  to the author (that's me).  I'd appreciate a reprint or e-copy
C	  of any paper where my code has been used.
C
C	- You can not sell it or use it in any way to get money out of it
C	  (save for the research grants you might get thanks to it)
C
C	- I'm making this code available to the best of my knowledge.
C	  It has been tested as throughly as possible. 
C            It work for Windows and Linux platforms.  However, it's
C	  your responsibility to check for its accuracy in the environment
C	  and size/shape ranges you use. 
C
C	- Disclaimer:  The usual stuff. I claim no responsibility for any
C	  misuse or damage arising from its use.  If this program erases your
C	  hard disk, makes your computer cry out for Linux, 
C	  makes you lose your boy/girlfriend, melts your TV set, spends out
C	  your lifesaving or forces you to get to decaf coffee, it's your
C	  problem.  For your tranquility, no such side effects has been
C	  reported ... yet.
C
C  WHILE THE COMPUTER PROGRAMS HAVE BEEN TESTED FOR A VARIETY OF CASES,
C  IT IS NOT INCONCEIVABLE THAT THEY MAY CONTAIN UNDETECTED ERRORS. ALSO,
C  INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
C  VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
C  THE AUTHORS DISCLAIM ALL LIABILITY FOR
C  ANY DAMAGES THAT MAY RESULT FROM THE USE OF THE PROGRAMS. 
C  I would highly appreciate informing me of any problems encountered 
C  with this code.  
C
C	  Improvements, comments, additions: all welcome at 
C
C	                   wavescattering@yahoo.com
C
C =====================================================================

      PROGRAM effe2p

C Variables declared but never referenced:
C    YNTEST declared at line 12 file effe.f

C Variables set but never used:
C    LMX set at line 20 file effe.f
C    TOL set at line 90 file effe.f
C Variable ZEPS is used before its value has been defined
C--------/---------/---------/---------/---------/---------/---------/--
C  PROGRAM TO CALCULATE EFFECTIVE REFRACTIVE INDEX USING THE 
C  RENORMALIZED MAXWELL-GARNETT FORMULA FOR A (COATED) SPHERE
C
C  		 f77 -g -C effe.f -o rneffe  
C--------/---------/---------/---------/---------/---------/---------/--
                                                                                   
      implicit none       
      integer NOUT,LMX,LCS,ILCS,NMAT,NFIN,NSTEP,NPC
      real*8 TOL
      logical ynperfcon,ynperfconv,yntest,ynbrug 
      COMPLEX*16 ZEPS0,CCEPS,CSEPS  
      COMPLEX*16 ZMU0,CCMU,CSMUS 
      character*1 ync,yncv

C ::: output unit number 
      PARAMETER (NOUT=60)
c number of spherical harmonics used
      PARAMETER (lmx=1)  
c number of different phases
      parameter (npc=2)          
*      
* If particle is coated, ync=y, otherwise ync=n (lower case only!)
      parameter (ync='n')
*
* ynperfcon=.true. if core is a perfect conductor, otherwise
* ynperfcon=.false.
*
      PARAMETER (ynperfcon=.false.)
*    
c number of coatings
      parameter (lcs=1)
c The coating layer to which material data are read in
      parameter (ilcs=1)
c if coated, the ratio 'core radius/particle radius'
c   (If lcs.ne.1, program is singular for rff=0. and 1.) - use homogeneous
c    particle instead!!!
c       PARAMETER (rff=0.95d0)
c background dielectric constant
      PARAMETER (ZEPS0=1.D0**2)
c
c background permeability
      PARAMETER (ZMU0=1.D0**2)
c
c Temporarily, before the coated part is finished,
c read in particle core dielectric function to CSEPS
c particle (core) dielectric constant (depending whether lcs=1 or lcs=2)
c If LCS.EQ.1, CCEPS is particle dielectric constant       
c      PARAMETER (CCEPS=(1.45D0,0.d0)**2)     !SiO2
	PARAMETER (CCEPS=(-50.84D0,0.762d0))    !JAP89_5774   ellipsoid for ld=633
c	PARAMETER (CCEPS=(-2.03D0,0.602d0))    !JAP89_5774   sphere for ld=354
C >>>     SPHERE (OUTER SHELL SCATTERER) PERMITTIVITY                  <<<
*  n(silica)=1.45  <--->    EPS(1)=2.1025D0
*  n(ZnS)=2.       <--->    EPS(1)=4.D0

      PARAMETER (CSEPS=(1.45d0,0.d0)**2)
c	PARAMETER (CSEPS=(-50.84D0,0.762d0))    !=(-10.84D0,0.762d0))    JAP89_5774

c If LCS.EQ.1, CCMU is particle permeability      

	PARAMETER (CCMU=(1.D0,0.d0))  

C >>>     SPHERE (OUTER SHELL SCATTERER) PERMEABILIITY                  <<<

	PARAMETER (CSMUS=(1.D0,0.d0)) 

* material code number 
c   NMAT=0             dispersionless dielectric                              
c   NMAT=1             Drude metal
c   NMAT=2             Ag
c   NMAT=3             Au
c   NMAT=4             ZnS 
c   NMAT=5             Cu
c   NMAT=6             Al
c   NMAT=7             Pt
c   NMAT=8             Si 
*
      PARAMETER(NMAT=1)   
*
c Temporarily option for reading of the real data for the dielectric constant 
c The number of the entries in a material data file to be read below
c          AGC.DAT                NFIN=73       ! from Palik 
c          Audat.dat              NFIN=66       ! from Palik
c          Au_2dat.dat            NFIN=76       ! from JAW  
c          Au*new.dat             NFIN=142
c          Cudat.dat              NFIN=47       ! from Palik
c          Aldat.dat              NFIN=80       ! from Palik
c          Ptdat.dat              NFIN=53       ! from Palik
c          Nidat.dat              NFIN=68       ! from Palik
c          Sidat.dat              NFIN=291
c          sieps.dat              NFIN=66
*
      PARAMETER (NFIN=66)
* 
C ::: relative error allowed for the TCS. If the convergence
*     within TOL is not reached, program issues warning
      PARAMETER (TOL=7.d-2)
*
c If ynbrug=.true., performs Bruggeman approximation for ZEPS1. Otherwise
c ynbrug=false.
      parameter (ynbrug=.false.)  
******************************************************************
c Declarations:

      integer ikl,ieps,istep,ipc
      REAL*8 rmuff,lambda,ff(npc),ffb,fft,pi,filfrac,enw,xstep,   !ffb,filfrac is for Bruggemann
     & RMF(lcs),rff(lcs),RMUF(npc),rsnm(npc)  
      real*8 omf(NFIN),omxf,reepsz,plasma,omxp
      real*8 delo,omega0,omega,omegal,p1,p2
      COMPLEX*16 ci,cone,xs,zeff,zneff,zmuff,z1,z2,zx,zy,
     & zalphae(npc),zalpham(npc),ca(npc),cb(npc),cc(npc),cd(npc),
     & cam(npc),cbm(npc),ccm(npc),cdm(npc)
	complex*16 ceps1(NFIN),ZEPS1,zeps(lcs+1),zmu(lcs+1)
      COMPLEX*16 TMT(4,3,3)
*
      COMMON /TOSPHERECR/ rmuff
      COMMON /TOSPHERECH/ yncv
      COMMON /TOSPHERECL/ ynperfconv
*
* From here to spheref
C
C--------/---------/---------/---------/---------/---------/---------/--

c Data:
      DATA PI/3.141592653589793d0/
      DATA ci/(0.d0,1.d0)/,cone/(1.d0,0.d0)/ 
*---------------------------------------------------------------
CCCCCCCCCCCCCCCCCC    Assignement of common variables CCCCCCCCCCC
*
      ynperfconv=ynperfcon
      yncv=ync
*****************************************************************

c      write(6,*)'The (coated) sphere radius in nanometers'
c      read(5,*) rsnm

	rsnm(1)=446d0    !almost close-packed LiTaO_{3} spheres 
        p1=0.5d0         !probability of sphere 1
	rsnm(2)=358d0/sqrt(zeps0)    !n-type Ge spheres 
        p2=1-p1        !probability of sphere 2

* sphere radius in the units, where
* the side ALPHA of a conventional unit cell of a cubic lattice
* is equal to 2 (then OMEGA=PI*ALPHA/LAMBDA). Various examples 
* of the sphere radius for different sphere filling fractions 
* and different lattices are listed below:


      RMUF(1)=rsnm(1)*2.D0/1264d0 
      ff(1)=p1*2.d0*pi*RMUF(1)**3/3.d0

      RMUF(2)=rsnm(2)*2.D0/1264d0 
      ff(2)=p2*2.d0*pi*RMUF(2)**3/3.d0
      fft=ff(1)+ff(2)              !total filling fraction of spheres


*oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
* Scanning over frequency interval:
*-------------------------------------------------
*
ct      write(6,*)'READ INITIAL (MAXIMAL) LAMBDA (in the vacuum) in nm'
ct      read(5,*) lambda
cv	lambda=633.d0                      !JAP89_5776 for gold ellipsod
c	lambda=354.d0                      !JAP89_5776 for silver sphere
cc      lambda=500d0

	lambda=13209d0
	       !299792.458d0/26.7d0/0.85d0  !lambda_T=11228,1819
	                                    !lambda_L= 6392,1632
*
*
* size parameter is customarily defined as the ratio of
* circumference of particle to the wavelength in the host medium
* in which the particle is embedded
*                    x=kr=\sg a=2*pi*r/lambda
*      xs=2.d0*pi*rs*dble(sqrt(zeps0))/lambda
* convert lambda to the lambda in vacuum:
c      lambda=lambda*dble(sqrt(zeps0))
c  omega=2.d0*pi*rs/(lambda*rmuf)=xs/(rmuf*dble(sqrt(zeps0))),
c where rs is the particle radius (in nm) and  lambda 
c is the wavelengths (in nm)
c in the vacuum:

         omega=2.d0*pi*rsnm(1)/(lambda*rmuf(1))
         omega0=omega

c      write(6,*)'Read omega'
c      read(5,*) omega
c      WRITE(6,*)'omega=', omega
c      if (1.eq.1) nstep=0             ! temporarily
c      delo=0.d0
c      omega0=omega
c      go to 11
* 
* Option for omega input:
c      write(6,*)'Read omega ='
c      read(5,*) omega
*
c      xs=RMUF*omega*dble(sqrt(zeps0))
*
c       write(6,*)'Equiv. size parameter x=2*pi*rs*n_0/lambda=',xs
*
ct      write(6,*)'Scan down to wavelength (in nm)'
ct      read(5,*) enw
c      enw=500
	enw=11224d0               !299792.458d0/46.9d0/0.95d0

ct      write(6,*)'Scanning step in nm'
ct      read(5,*) xstep

      xstep=2
*
C ::: number of steps on frequency interval:
      if (lambda.eq.enw) then
       NSTEP=0
       delo=0.d0
      else
       NSTEP=(lambda-enw)/xstep
C ::: width of the searched frequency interval
       ENW=2.d0*pi*rsnm(1)/(enw*rmuf(1))
       enw=enw-omega0
       delo=enw/dble(nstep)
      end if
*

      ZEPS(1)=cceps
      if((lcs.gt.1).and.
     & ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.lcs))))
     & zeps(lcs)=cseps
      zeps(lcs+1)=zeps0

	ZMU(1)=ccmu
      if((lcs.gt.1).and.
     & ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.lcs))))
     & zmu(lcs)=csmus
      zmu(lcs+1)=zmu0

*                  --------------------------------
* output initial statements


      OPEN(UNIT=NOUT,FILE='effmed.dat')
      rewind(NOUT)
      WRITE(NOUT,*)'#EFFECTIVE MEDIUM PARAMETERS'
      write(nout,*)
      write(nout,*)'The effective dielectric constant obtained in the 
     & non-quasi-static case is complex even in the case of 
     & non-absorbing spheres embedded in a non-absorbing medium. 
     & Therefore, except for the quasi-static case, the imaginary 
     & of the effective dielectric constant of the extended 
     & MG theories cannot be used to evaluate the absorption 
     & or the heating rate.'
      write(nout,*)
	write(nout,*)'#Spheres filling fraction=', fft
	write(nout,*)'#Sphere 1 radius in nm=', rsnm(1)
	write(nout,*)'#Sphere 1 radius rmuf in relative units=', rmuf(1)
        write(nout,*)'#Sphere 2 radius in nm=', rsnm(2)
        write(nout,*)'#Sphere 2 radius rmuf in relative units=', rmuf(2)
      if (.not.ynperfcon) then
        write(nout,*)'#Material number =',NMAT
      else 
        write(nout,*)'#Perfectly conducting sphere (core)'
      end if   
      if (ynbrug) write(nout,*)'#Bruggeman approximation performed'
      if (ync .eq.'n') then
        write(nout,*)'#Homogeneous particle'
      else if (ync.eq.'y') then
       write(nout,*)'#Coated particle with (including core)',
     & LCS,'shells'
      end if 
      if (nmat.ge.1) write(nout,*)'#Dispersive layer number=', ilcs
      write(nout,*) 
      write(nout,*)'#host dielectric constant=', zeps(lcs+1)
      if ((lcs.eq.1).and.(nmat.eq.0)) 
     &     write(nout,*)'#sphere diel. const.=', zeps(1)  
      if (lcs.ge.2) then
      write(nout,*)'#core radius/sphere radius =',rff(1)
      if ((ilcs.ne.1).or.(nmat.eq.0))  
     &  write(nout,*)'#sphere core diel. const.=', zeps(1)  
      if ((ilcs.ne.lcs).or.(nmat.eq.0))  
     &  write(nout,*)'#coating diel. const.=',zeps(lcs)
C--------/---------/---------/---------/---------/---------/---------/-- 
      end if  
      WRITE(NOUT,*)
      WRITE(NOUT,*)'#In columns: lambda_0 and the respective Re and Im'
      WRITE(NOUT,*)'#parts of eff. diel. constant, eff. mu, and n_eff'
      write(nout,*)
C--------/---------/---------/---------/---------/---------/---------/-- 



C--------/---------/---------/---------/---------/---------/---------/--
* Checking set up:
*
      if (NMAT.GT.1) write(6,*)'Real material data 
     & are to be provided'
*
      if ((ync.eq.'y'.and.lcs.eq.1).or.(ync.eq.'n'.and.lcs.ne.1)) then
      write(6,*)'Check compatibility of YNC and LCS'
      stop
      end if
      
C--------------------------------------------------------------------

*
      if (lcs.ge.2) then                  !coated particle
*
cc      write(6,*)'Read equal-volume-core radius in nm'
cc      read(5,*) rff(1)
cc      rff(1)=204.9d0
cc      rff(1)=rff(1)/rs

      write(6,*)'Coated sphere shell radii r(l) labelled from 1 for the 
     & inner core till LCS for the outer shell
     &   ===> r(lcs) is the sphere radius'
C--------/---------/---------/---------/---------/---------/---------/--
      do ikl=1,lcs-1
*
      write(6,*)'Read r(l) for l=',ikl
      read (5,*) rff(ikl)
      rff(ikl)=rff(ikl)/rsnm(ipc)

c      rff(1)=0.75d0
      rmf(ikl)=rff(ikl)*rmuff

      IF ((IKL.GE.2).and.(IKL.LT.LCS)) THEN

      if ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.ikl))) then
      write(6,*)'Read the lth-sphere layer diel. const. for l=',ikl
      read (5,*) zeps(ikl)
      end if
*
      write(6,*)'Read the lth-sphere layer magnetic permeability
     &   for l=',ikl
      read (5,*) zmu(ikl)

      END IF

	enddo                 !IKL.GE.2
*
      end if                !lcs.ge.2



*****************************   ZEPS1  ***********************************
*                  --------------------------------       
* Sphere optical constants in the case of a dispersion
* READING IN MATERIAL  DATA:
* Reading real silver data according to Palik's  book
* requires reading data files OMF and CEPS1 of dimension NFIN
* OMF is reepsz/omega and CEPS1 contains the sphere EPS
*                       material constant reading:
*
      if (nmat.eq.2) then            ! silver data
      OPEN(UNIT=30,FILE='agc.dat')           
      rewind(30)
        do ieps=1,nfin
          read(30,*) omf(ieps),ceps1(ieps)
        enddo
       close(30)
       
      else if (nmat.eq.3) then        ! Gold data 

c      OPEN(UNIT=30,FILE='Au293Knew.dat')       !Gold data for different T
      OPEN(UNIT=30,FILE='Audat.dat')          !Gold data in nm
      write(6,*)'Gold particles'
      rewind(30)
        do ieps=1, nfin
          read(30,*) omf(ieps),ceps1(ieps)
c          omf(ieps)=2.d0*pi*rsnm*omf(ieps)/(1240.d0*rmuf)
          omf(ieps)=2.d0*pi*rsnm(1)/(omf(ieps)*rmuf(1))
        enddo
       close(30)

cc      else if (nmat.eq.4) then          
      
      else if (nmat.eq.5) then        ! Copper data

      OPEN(UNIT=30,FILE='Cudat.dat')          !Copper data in nm
      write(6,*)'Copper particles'
      rewind(30)      
        do ieps=1, nfin
          read(30,*) omf(ieps),ceps1(ieps)
          omf(ieps)=2.d0*pi*rsnm(1)/(omf(ieps)*rmuf(1))
        enddo      
      close(30)     

      else if (nmat.eq.6) then        ! Aluminium data 

      OPEN(UNIT=30,FILE='Aldat.dat')          !Aluminium data in nm
      write(6,*)'Aluminum particles'
      rewind(30)      
        do ieps=1, nfin
          read(30,*) omf(ieps),ceps1(ieps)
          omf(ieps)=2.d0*pi*rsnm(1)/(omf(ieps)*rmuf(1))
        enddo      
      close(30) 

      else if (nmat.eq.7) then        ! Platinum data 

      OPEN(UNIT=30,FILE='Ptdat.dat')          !Platinum data in nm
      write(6,*)'Platinum particles'
      rewind(30)      
        do ieps=1, nfin
          read(30,*) omf(ieps),ceps1(ieps)
          omf(ieps)=2.d0*pi*rsnm(1)/(omf(ieps)*rmuf(1))
        enddo      
      close(30) 
            
      end if                      ! material constant reading      

*********************
*                     --------------------------------
* begin main scanning loop:

      do 200 istep=1,nstep+1

      omega=omega0 + dble(istep-1)*delo
      

C--------/---------/---------/---------/---------/---------/---------/-- 

      do 100 ipc=1,npc

	xs=RMUF(ipc)*omega*dble(sqrt(zeps0))  !Equiv. size parameter

      rmuff=RMUF(ipc)
      rmf(lcs)=rmuff

      if ((nmat.eq.0).or.(ynperfcon)) go to 8      !dispersionless dielectric 
                                                   !or ideal metal
* In case of a dispersion, ZEPS1 is modified.
* For ideal Drude metal
*     plasma=2.d0*pi*particle radius in nm/(lambda_z in nm*rmuf)
* where lambda_z is the wavelength for which Re eps_s=0.

C      reepsz=2.d0*pi*rsnm(1)/(324.269d0*rmuf(1))   !setting to Ag plasma frequency
C	plasma=reepsz

* omega_p=omega_T of 2*pi*26.7 THZ corresponds to
* in local frequency units



      IF (NMAT.EQ.1) THEN              !Material decision IF - Drude metal

* c = 2.99792458d8 m/s
* THz is equivalent to $\ld=2.99792458\times 10^{-4} m$, or to $299,792458$ $\mu$m.

	 plasma=2.d0*pi*rsnm(1)*26.7d0/(299792.458d0*rmuf(1))    !omega_T
	 omegal=2.d0*pi*rsnm(1)*46.9d0/(299792.458d0*rmuf(1))    !omega_L

	if (ipc.eq.1) then

c	 omegal=2.d0*pi*rsnm(1)*46.9d0/(299792.458d0*rmuf(1))

      zeps1=13.4d0*(1.d0 + (omegal**2-plasma**2)/(plasma**2-omega**2))

	else if (ipc.eq.2) then

c        plasma=omegal     
        omxp=plasma*1.d0/omega                !=plasma/omega

        zeps1=15.8d0*(1.d0 - omxp**2)
	
	end if
      
	go to 5
*
      ELSE IF (nmat.eq.4) then             !Material decision IF - ZnS
*
       filfrac=0.62d0         ! filfrac of ZnS in ZnS core
       call  znsrefind(LAMBDA,FILFRAC,zeps1) 
       go to 5 
*
      ELSE IF (NMAT.EQ.2) THEN         !Material decision IF - Ag

c >>> real material data:           !silver 
*                         lambda_z=323.83d0
*                         lambda_p=164.d0
* When real material data are used, 
* reepsz differs from plasma!!! The plasma wavelength is 
* calculated below: 

       plasma=reepsz*7.2d0/3.8291d0   !setting to Au plasma frequency

* security trap - remainder (not optimized!)
      omxf=omega/reepsz
      if (omxf.gt.omf(1)) then
       write(6,*)'Calculation of has to stop with'
       write(6,*)' OMF(1)'
       write(6,*)' OMXF=', omxf
       stop
      end if

      if (omxf.lt.omf(nfin)) then
        omxp=plasma/omega
        zeps1=1.d0-omxp**2/(1.d0+ci*plasma/(144.d0*omega))
* damping coefficient for silver is plasma/144 where plasma is different from
* the Re eps zero crossing at 3.8291 eV according to Palik!!!
       go to 5
      else if (omxf.eq.omf(1)) then
       zeps1=ceps1(1)
       go to 5
      else
      do ieps=2,nfin
* data file ordered with the increased wavelength
* omxf increases in the loop and is oriented opposite to the data file
       if (omxf.gt.omf(ieps)) then     ! linear interpolation
       zeps1=ceps1(ieps)+(omxf-omf(ieps))*(ceps1(ieps-1)-ceps1(ieps))
     1 /(omf(ieps-1)-omf(ieps))
       go to 5
       end if 
      enddo
       end if 
*
*
      ELSE IF ((NMAT.GE.3).or.((nmat.ge.5).and.(nmat.le.7))) then   !Material decision IF 
                                                                    !Au,Cu,Al,Pt 
c >>> 
* data file ordered with the decreased wavelength
* omega increases in the loop and is oriented along the data file
*
      if ( (omega.lt.omf(1)).or.(omega.gt.omf(nfin)) ) then
cc       write(6,*)'Material data not available for this wavelength'
cc       stop
*
      call sordalc(NMAT,lambda,ZEPS1)
      go to 5
*
      end if
*
      if (omega.eq.omf(nfin)) then
       zeps1=ceps1(nfin)
       go to 5
      else 
      do ieps=1,nfin-1
       if (omega.lt.omf(ieps+1)) then     ! linear interpolation
       zeps1=ceps1(ieps)+(omega-omf(ieps))*(ceps1(ieps+1)-ceps1(ieps))
     1 /(omf(ieps+1)-omf(ieps))
       go to 5
       end if 
      enddo
      end if

       END IF                  ! END of Material decision IF 

*The end of reading real data according to Palik's  book   
*_____________________________________
* activate Bruggeman:

  5   if (ynbrug) then
      ffb=0.8d0 
      z1 = (3.d0*ffb-1.d0)*zeps1+(2.d0 - 3.d0*ffb)*zeps0
      z2 =  sqrt(z1*z1 + 8.d0*zeps1*zeps0)
*
       if (IMAG(z2).GE.0.0) then
         zeps1= (z1 + z2)/4.d0
       else            
         zeps1= (z1 - z2)/4.d0
       end if
       end if

      zeps(ilcs)=zeps1
*______________________________________

  8   lambda=2.d0*pi*rsnm(1)/(omega*RMUF(1))
*
      call spheref(1,lcs,lambda,rsnm(ipc),rmuf(ipc),zeps,zmu,tmt)  
C--------/---------/---------/---------/---------/---------/---------/--
C >>> lmax,lcs,lambda,rsnm(ipc),rmf(ipc),zeps,zmu
C <<< tmt
C ===============
C  This routine calculates the single sphere scattering properties
C          (including coated spheres)
C
C Outputs the elements of the diagonal T matrix
C
C                         |  TMT(1,*) |       0     |
C                 TMT  =  | ----------+-------------|
C                         |     0     |  TMT(2,*)   |
C
C    TMT(1,*) terms corresponds to TEE scattering matrices
C    TMT(2,*) terms corresponds to TMM scattering matrices
C
C    TMT's equal to i*sin(eta)*exp(i*eta), where eta is a phase-shift
C--------/---------/---------/---------/---------/---------/---------/--

c dielectric polarization factor:

      zalphae(ipc)=-ci*3.d0*tmt(1,1,1)/(2.d0*xs**3)

      ca(ipc)= zeps0*(cone +2.d0*zalphae(ipc)*fft)
      cb(ipc)= cone - zalphae(ipc)*fft
      cc(ipc)= zeps0*(cone + 2.d0*zalphae(ipc)*fft)
      cd(ipc)= 2.d0*cone -2.d0*zalphae(ipc)*fft


c magnetic polarization factor:

      zalpham(ipc)=-ci*3.d0*tmt(2,1,1)/(2.d0*xs**3)

	zmuff=zmu0*(cone + 2.d0*zalpham(ipc)*ff(ipc))/
     & (cone - zalpham(ipc)*ff(ipc))
	write(6,*)'MG zmuff=', zmuff

      cam(ipc)= zmu0*(cone +2.d0*zalpham(ipc)*fft)
      cbm(ipc)= cone - zalpham(ipc)*fft
      ccm(ipc)= zmu0*(cone + 2.d0*zalpham(ipc)*fft)
      cdm(ipc)= 2.d0*cone -2.d0*zalpham(ipc)*fft


 100  continue              ! ipc-loop over different spheres

c Luo two-phase Maxwell-Garnett formula

      zx=-(p1*(ca(1)*cd(2)-cb(1)*cc(2))+p2*(ca(2)*cd(1)-cb(2)*cc(1)))/
     & (p1*cb(1)*cd(2)+p2*cb(2)*cd(1))
      zy=-(p1*ca(1)*cc(2)+p2*ca(2)*cc(1))/
     & (p1*cb(1)*cd(2)+p2*cb(2)*cd(1))
C--------/---------/---------/---------/---------/---------/---------/--


      zeff=(-zx + sqrt(zx**2-4.d0*zy))/2.d0
      if (imag(zeff).lt.0) zeff=(-zx - sqrt(zx**2-4.d0*zy))/2.d0


      zx=-(p1*(cam(1)*cdm(2)-cbm(1)*ccm(2))+p2*(cam(2)*cdm(1)
     & -cbm(2)*ccm(1)))/(p1*cbm(1)*cdm(2)+p2*cbm(2)*cdm(1))
      zy=-(p1*cam(1)*ccm(2)+p2*cam(2)*ccm(1))/
     & (p1*cbm(1)*cdm(2)+p2*cbm(2)*cdm(1))
C--------/---------/---------/---------/---------/---------/---------/--


      zmuff=(-zx + sqrt(zx**2-4.d0*zy))/2.d0
      if (imag(zmuff).lt.0) zmuff=(-zx - sqrt(zx**2-4.d0*zy))/2.d0

 
	zneff=sqrt(zeff*zmuff)
	if (imag(zneff).lt.0)	zneff=-sqrt(zeff*zmuff)

C
	write(6,*)
      write(6,*)
	write(6,*)'LAMBDA=',lambda
      write(6,*)
	write(6,*)'QUASISTATIC VALUES:'
      write(6,*)'effective dielectric constant e_eff=', zeff
	write(6,*)'effective permeability mu_eff=', zmuff
      write(6,*)'effective complex refractive index n_eff=',zneff
C
	write(nout,1105) lambda, zeff, zmuff, zneff

C

 200  continue

*
      close(nout)

      write(6,*)
      write(6,*)'Effective medium parameters versus wavelength in
     & effmed.dat'
      write(6,*)
      write(6,*)'The effective dielectric constant obtained in the 
     & non-quasi-static case is complex even in the case of 
     & non-absorbing spheres embedded in a non-absorbing medium. 
     & Therefore, except for the quasi-static case, the imaginary 
     & of the effective dielectric constant of the extended 
     & MG theories cannot be used to evaluate the absorption 
     & or the heating rate.'
C--------/---------/---------/---------/---------/---------/---------/--
 
 1105 FORMAT(F12.3,4X,F10.4,2X,F10.4,4X,F10.4,2X,F8.5,4X,F10.4,2X,F10.4)

      END
C********************************************************************


      subroutine spheref(lmax,lcs,lambda,rsnm,rmf,zeps,zmu,tmt)  
C--------/---------/---------/---------/---------/---------/---------/--
C >>> lmax,lcs,lambda,rsnm,rmf,zeps,zmu
C <<< tmt
C
C  Outputs the elements of the diagonal T matrix
C
C                         |  TMT(1,*) |       0     |
C                 TMT  =  | ----------+-------------|
C                         |     0     |  TMT(2,*)   |
C
C    TMT(1,*) terms corresponds to TEE scattering matrices
C    TMT(2,*) terms corresponds to TMM scattering matrices
C
C    TMT's equal to i*sin(eta)*exp(i*eta), where eta is a phase-shift
C ===============
C  This routine calculates the single sphere scattering properties
C          (including coated spheres)
C
C
C k_l length in units (2*PI/A=) PI:    xkl= 0.8660254037844386d0
C
C KMT ... elements of the K matrix
C ALPHA, BETA ... arrays of the electric and magnetic phase shifts
C 
C Partial wave expansion is used, which is badly convergent 
C for large size parameters $x> 100$. In numerical applications, the 
C series is to be cut off after 
C          LMAX (LMX parameter here) \approx x+4x^{1/3}+2$.
C In the case of LMX=50 this means that x <= 35.
C If one wants to observe ripples, the cutoff for a given x has to 
C be even larger
C---------------------------------------------------------------------                                                                                     
      implicit none       
      integer LMX,LCS,LMAX,LMAXD1,LMTD
      real*8 TOL
      character*1 ync
      logical ynperfcon

c Parameters:
c number of spherical harmonics used
      PARAMETER (lmx=1,LMAXD1=LMX+1,LMTD=LMAXD1*LMAXD1-1)

c maximal number of sphere coatings
cc      parameter (lcs=3)
c
C ::: relative error allowed for the TCS. If the convergence
*     within TOL is not reached, program issues warning
      PARAMETER (TOL=1.d-6)
*
C--------/---------/---------/---------/---------/---------/---------/--
c Declarations:
      integer i,j,l,ikl,ij,ij1,m
      real*8 ZEPS0
      real*8 rmf(lcs),rmuf,xs,rsnm,pi
      real*8 omega,lambda
*
      COMPLEX*16 ci,czero,zeps(lcs+1),zmu(lcs+1),cqeps(2),cqmu(2)
      COMPLEX*16 RX(2),SG(2),zs,zz1,zz2
      COMPLEX*16 KMT(2,lmx),TMT(4,LMTD,LMTD),alpha(lmx),beta(lmx)
*
      COMPLEX*16 cm(lcs,lmx),dm(lcs,lmx),ce(lcs,lmx),de(lcs,lmx)
      COMPLEX*16 tt1(2,2,lmx,2),tt2(2,2,lmx,2)
      COMPLEX*16 AM(lmx),AE(lmx),BM(lmx),BE(lmx)
*
      COMPLEX*16 JL(0:lmx),NL(0:lmx)
      COMPLEX*16 DRJL(0:lmx),DRNL(0:lmx)
      COMPLEX*16 UL(2,0:lmx),VL(2,0:lmx)
      COMPLEX*16 DRUL(2,0:lmx),DRVL(2,0:lmx)
      COMPLEX*16 ZARTAN   
* 
      external ZARTAN
*
*
      COMMON /TOSPHERECR/ rmuf                  !rmuf(ipc)
      COMMON /TOSPHERECH/ ync
      COMMON /TOSPHERECL/ ynperfcon
*
* From the main here
*---------------------------------------------------------------
c Data:
      DATA PI/3.141592653589793d0/
      DATA ci/(0.d0,1.d0)/,czero/(0.D0,0.D0)/
C--------/---------/---------/---------/---------/---------/---------/--
* Initialization:

      do ij=1,lmtd
       do ikl=1,lmtd
         tmt(3,ikl,ij)=czero
         tmt(4,ikl,ij)=czero

           if (ikl.ne.ij) then
             tmt(1,ikl,ij)=czero
             tmt(2,ikl,ij)=czero
           end if

       enddo
      enddo

C--------/---------/---------/---------/---------/---------/---------/--
* Checking set up:

      if ((ync.eq.'y'.and.lcs.eq.1).or.(ync.eq.'n'.and.lcs.ne.1)) then
      write(6,*)'Check compatibility of YNC and LCS'
      stop
      end if

C--------------------------------------------------------------------
* Reading in the input data: 
c close packed :
c       RMUF=1.d0/DSQRT(2.D0)
*
* size parameter is customarily defined as the ratio of
* circumference of sphere to the wavelength in the host medium
* in which the sphere is embedded
*                    x=kr=\sg a=2*pi*r/lambda
*      xs=2.d0*pi*rs*dble(sqrt(zeps0))/lambda
* convert lambda to the lambda in vacuum:
c      lambda=lambda*dble(sqrt(zeps0))
c  omega=2.d0*pi*rsnm/(lambda*rmuf)=xs/(rmuf*dble(sqrt(zeps0))),
c where rs is the sphere radius (in nm) and  lambda is the wavelengths (in nm)
c in the vacuum:
* 
* Option for omega input:
c      write(6,*)'Read omega ='
c      read(5,*) omega
*
         zeps0=zeps(lcs+1)
         omega=2.d0*pi*rsnm/(lambda*rmuf)    !=2.d0*pi*rsnm(ipc)/(lambda*rmf(ipc))

c      xs=2.d0*pi*rsnm*sqrt(zeps0)/lambda 
c      xs=RMUF*omega*dble(sqrt(zeps0))

*
c       write(6,*)'Size parameter x=2*pi*rs*n_0/lambda=', xs
*
*                  --------------------------------
*______________________________________
*

      if (.not.ynperfcon) then
      
      ij1=1

      do l=1,lmx
      AM(l)=dcmplx(1.d0,0.d0)
      AE(l)=dcmplx(1.d0,0.d0)
      BM(l)=dcmplx(0.d0,0.d0)
      BE(l)=dcmplx(0.d0,0.d0)
      enddo
           
      else if (ynperfcon) then
            
      CQEPS(2)=SQRT(ZEPS(2))
      CQMU(2)=SQRT(ZMU(2))
      SG(2)=omega*SQRT(ZEPS(2)*ZMU(2))
      RX(1)=SG(2)*RMF(1)
*
*
      call gnzbess(RX(1),LMX,jl,drjl,nl,drnl)
*      
      DO 10 L=1,lmx
C >>> (AS 10.1.22):
      UL(1,L)=RMF(1)*JL(L)
      VL(1,L)=RMF(1)*NL(L)
      DRJL(L)=SG(2)*DRJL(L)
      DRNL(L)=SG(2)*DRNL(L)
      DRUL(1,L)=JL(L)+RMF(1)*DRJL(L)
      DRVL(1,L)=NL(L)+RMF(1)*DRNL(L)
      AM(l)= NL(L)                ! cm(1,l)
      BM(l)=-JL(L)                ! dm(1,l)
      AE(l)= DRVL(1,L)            ! ce(1,l)
      BE(l)=-DRUL(1,L)            ! de(1,l)

* cf. Jackson 1962, p. 571, Eqs. (16.147); 
*                        B/A should yield -tan(phase shift)


  10  continue 
 
      if (lcs.eq.1) go to 30
      
      ij1=2 
      
      end if             
C********************************************************************
c Execution:
* Calculation of the phase shifts
*
      DO 28 j=ij1,lcs

      CQEPS(1)=SQRT(ZEPS(j))
      CQMU(1)=SQRT(ZMU(j))

      SG(1)=omega*SQRT(ZEPS(j)*ZMU(j))
*
      CQEPS(2)=SQRT(ZEPS(j+1))
      CQMU(2)=SQRT(ZMU(j+1))
      SG(2)=omega*SQRT(ZEPS(j+1)*ZMU(j+1))
*
      DO 25 I=1,2
*
      RX(I)=SG(I)*RMF(j)
c      WRITE(6,*)'i, rx(i)=', i, rx(i)
C >>>
*
      call gnzbess(RX(I),lmx,jl,drjl,nl,drnl)
*
c      write(6,*)'jl=', jl 
      DO 15 L=1,lmx
C >>> (AS 10.1.22):
      UL(I,L)=RMF(j)*JL(L)
      VL(I,L)=RMF(j)*NL(L)
      DRJL(L)=SG(I)*DRJL(L)
      DRNL(L)=SG(I)*DRNL(L)
      DRUL(I,L)=JL(L)+RMF(j)*DRJL(L)
      DRVL(I,L)=NL(L)+RMF(j)*DRNL(L)

  15  continue

  25  CONTINUE
*
c      write(6,*)'ul=', ul 
*
C >>>  END OF THE LOOP TO ASSIGN VALUES OF BESSEL FUNCTIONS
C      JL and NL start to oscillate after RX.GT. approx 2.5
C********************************************************************
C           Transfer matrix for a layered (coated) sphere
C********************************************************************
*
      do l=1,lmx
*
*   magnetic part
*
      tt1(1,1,l,1)= cqmu(1)*UL(1,L)
      tt1(1,2,l,1)= cqmu(1)*VL(1,L)
      tt1(2,1,l,1)= DRUL(1,L)/cqmu(1)
      tt1(2,2,l,1)= DRVL(1,L)/cqmu(1)
*
      tt2(1,1,l,1)= sg(2)*DRVL(2,L)/cqmu(2)
      tt2(1,2,l,1)= - sg(2)*cqmu(2)*VL(2,L)
      tt2(2,1,l,1)= - sg(2)*DRUL(2,L)/cqmu(2)
      tt2(2,2,l,1)= sg(2)*cqmu(2)*UL(2,L)
*
*   electric part
*
      tt1(1,1,l,2)=cqeps(1)*UL(1,L)
      tt1(1,2,l,2)=cqeps(1)*VL(1,L)
      tt1(2,1,l,2)=DRUL(1,L)/cqeps(1)
      tt1(2,2,l,2)= DRVL(1,L)/cqeps(1)
*
      tt2(1,1,l,2)= sg(2)*DRVL(2,L)/cqeps(2)
      tt2(1,2,l,2)= -sg(2)*cqeps(2)*VL(2,L)
      tt2(2,1,l,2)= -sg(2)*DRUL(2,L)/cqeps(2)
      tt2(2,2,l,2)= sg(2)*cqeps(2)*UL(2,L)
*
* m-part
*
      cm(j,l)=AM(l)*(tt2(1,1,l,1)*tt1(1,1,l,1)
     1 +tt2(1,2,l,1)*tt1(2,1,l,1))+BM(l)*(
     2 tt2(1,1,l,1)*tt1(1,2,l,1)+tt2(1,2,l,1)*tt1(2,2,l,1))
*
      dm(j,l)=AM(l)*(tt2(2,1,l,1)*tt1(1,1,l,1)
     1 +tt2(2,2,l,1)*tt1(2,1,l,1))+BM(l)*(
     2 tt2(2,1,l,1)*tt1(1,2,l,1)+tt2(2,2,l,1)*tt1(2,2,l,1))
*
* e-part
*
      ce(j,l)=AE(l)*(tt2(1,1,l,2)*tt1(1,1,l,2)
     1 +tt2(1,2,l,2)*tt1(2,1,l,2))+BE(l)*(
     2 tt2(1,1,l,2)*tt1(1,2,l,2)+tt2(1,2,l,2)*tt1(2,2,l,2))
*
      de(j,l)=AE(l)*(tt2(2,1,l,2)*tt1(1,1,l,2)
     1 +tt2(2,2,l,2)*tt1(2,1,l,2))+BE(l)*(
     2 tt2(2,1,l,2)*tt1(1,2,l,2)+tt2(2,2,l,2)*tt1(2,2,l,2))
*
      AM(l)=cm(j,l)
      BM(l)=dm(j,l)
      AE(l)=ce(j,l)
      BE(l)=de(j,l)
c      write(6,*) AM(l), BM(l)
c      write(6,*) AE(l), BE(l)
*
      enddo
*
  28  CONTINUE
c      write(6,*)'am=', am
c      write(6,*)'bm=', bm
C--------/---------/---------/---------/---------/---------/---------/--
C     ASSIGNING VALUES TO ELEMENTS OF THE K-MATRIX
C >>>
  30  CONTINUE
*
      DO 40 L=1,lmx

* In the following, one needs only phase shifts, so that division
* by SG(2) is omitted:
*
      KMT(1,L)=bm(l)/am(l)
      KMT(2,L)=be(l)/ae(l)
*
 40   CONTINUE
c      write(6,*)'kmt=', kmt
C********************************************************************

      DO 60 j=1,lmx

* >>> Extracting phase-shifts from KMT
*  >>> magnetic part:

         zs=-kmt(1,j)
         alpha(j)=zartan(zs)
* 
* Under normal circumstances, Im. part of a phase shift \geq 0 !!!)
*
      if(dimag(alpha(j)).lt.0.d0) then
      write(6,*)'dimag(alpha(j))=',dimag(alpha(j)) ,' is negative'
      write(6,*)'omega, j=', omega, j
      pause
      end if
c      write(6,*)'j, alpha(j)=', j, alpha(j)
*>>> electric part:

         zs=-kmt(2,j)
         beta(j)=zartan(zs)
* 
* Under normal circumstances, Im. part of a phase shift \geq 0 !!!)
*
      if(dimag(beta(j)).lt.0.d0) then
      write(6,*)'dimag(alpha(j))=',dimag(beta(j)) ,' is negative'
      write(6,*)'omega, j=', omega, j
      pause
      end if
c      write(6,*)'beta(j)=', beta(j)

  60  CONTINUE

      do 100 l=1,lmax

      zz1 = ci*sin(beta(l))*exp(ci*beta(l))
      zz2 = ci*sin(alpha(l))*exp(ci*alpha(l))

      do 100 m=-l,l
    
        ij=l*(l+1)+m
         
        tmt(1,ij,ij)= zz1 
        tmt(2,ij,ij)= zz2

 100  continue
*
      return
      end

C (C) Copr. 11/2003  Alexander Moroz
C********************************************************************

      function zartan(zs)
*--------/---------/---------/---------/---------/---------/---------/--
* For a complex argument using that 
*              atan (z)= -ln((1+iz)/(1-iz))*i/2.d0
* Hower, a direct application of this formula often
* results in a nonzero imaginary part of  $\arctan z$ even
* for a  purely real $z=x$.
* Therefore, in order to avoid this, in a numerical
* implementation, the complex logarithm is written
* explicitly as
*
*  ln\left(\fr{1+iz}{1-iz}\right)= (ar1,ar2)= log(
*  (1.d0 + abs(z)**2 -2.d0*imag(z))/1.d0+abs(z)**2 +2.d0*imag(z)))/2.d0
* +
* ci*atan(2.d0*dble(z)/(1-abs(z)**2).
*
* For a real z=x, |x|<1,  log here is purely imaginary and
* equals to i*atan(2x/(1-x**2))\equiv 2*i*atan x
*--------/---------/---------/---------/---------/---------/---------/--
      complex*16 zartan,zs
      real*8 xxs,xas,ar1,ar2,pi
      DATA PI/3.141592653589793d0/

      if (dimag(zs).eq.0.) then
        zartan=dcmplx(atan(dble(zs)),0.d0)
      return
      end if
         xas=abs(zs)**2
      ar1=log((1.d0+xas-2.d0*dimag(zs))/(1.d0+xas+2.d0*dimag(zs)))/2.d0
         xxs=dble(zs)
* special case:
         if(xas.eq.1.) then 
           if(xxs.ge.0.) then
            ar2=pi/2.d0
           else if (xxs.lt.0.) then
            ar2=-pi/2.d0
           end if
         zartan =dcmplx(ar2,- ar1)/2.d0
         end if

* remaining cases: 
         ar2=2.d0*xxs/(1.d0-xas)

         if(xas.lt.1.d0)  then     ! 1st and 4th quadrant
         zartan=dcmplx(atan(ar2),- ar1)/2.d0
         else if (xas.gt.1. .and. xxs.ge.0.) then       ! 2nd quadrant
         zartan=dcmplx(pi+atan(ar2),- ar1)/2.d0
         else if(xas.gt.1. .and. xxs.lt.0.) then        ! 3rd quadrant
         zartan=dcmplx(-pi+atan(ar2),- ar1)/2.d0
         end if
       return
       end
C (C) Copr. 1/1999  Alexander Moroz
C********************************************************************

     subroutine sordalc(NMAT,LAMBDA,ZEPS)
C----------------------------------------------------------------------
C        SUBROUTINE TO CALCULATE THE DIELECTRIC CONSTANT OF METALS
C             ACCORDING TO AN ARTICLE BY M. A. Ordal et al,
C  Optical properties of fourteen metals in the infrared and far infrared:
C   Al, Co, Cu, Au, Fe, Pb, Mo, Ni, Pd, Pt, Ag, Ti, V, and W,
C   Appl. Opt. {\bf 22}, 1099 (1983); ibid. {\bf 24}, 4493 (1985)
C  
C             f77 -g -check_bounds ordalc.f -o rnordalc
C          
C   omega=1/\lambda [cm^{-1}]  in Ordal [spectroscopic convention]
C
C   Re eps = - \fr{\om_p^2}{\om^2+\om_\tau^2}
C   Im eps =   \fr{\om_p^2 \om_\tau}{\om^3+\om \om_\tau^2}
C   1eV=1.24 \mu m
C----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER NMAT
      REAL*8 lambda,plasma,tau
      COMPLEX*16 zeps
C                       -------------------------------
      REAL*8 pi,x,y,omega
      DATA PI/3.141592653589793d0/
C   ---------      
C ::: speed of light in vacuum in nm/s
C      PARAMETER (c0=2.9927925d17)
C
C According to Table I of Ordal et al, Appl. Opt. {\bf 24}, 4493 (1985):
C                 (plasma=omega_plasma and tau=omega_tau below)
C
C          plasma[THz]/plasma[cm-1]  tau[THz]/tau[cm-1]   
C   Al            3570/119000            19.4/660          
C   Cu            1914/59600             8.34/73.2           
C   Au            2175/72800              6.5/215            
C   Ag            2175/72700             4.35/145        
C   Pt                /41500                 /558       
C
C ::: conversion factor between normal angular frequency and eV:
c      PARAMETER(XCN=4.13566727d-15)    

      if (nmat.eq.3) then           !Au
C ::: Convert the plasma frequency from Table I of Ordal et al
C ::: from [cm-1] to [nm-1] ===> conversion factor 10^{-7}
         PLASMA=72800.d-7            !d12 in Hz/d-7 in [nm-1]
C ::: Convert the tau frequency from Table I of Ordal et al
C ::: from [cm-1] to [nm-1] ===> conversion factor 10^{-7}
          TAU=215.d-7     
      else if (nmat.eq.5) then      !Cu
        PLASMA=59600d-7
        TAU=73.2d-7
      else if (nmat.eq.6) then      !Al
        PLASMA=119000d-7
        TAU=660.d-7
      else if (nmat.eq.7) then      !Pt 
        PLASMA=41150d-7
        TAU=558.d-7
      end if
C                       -------------------------------

c      write(6,*)'Read in wavelength in nm'
c      read(5,*) lambda
       omega=1.d0/lambda
*
       X=-plasma**2/(omega**2+tau**2)
       Y=tau*plasma**2/(omega**3+omega*tau**2)       
*
       zeps=dcmplx(X,Y)  
*
       END
*
C (C) Copr. 04/2002  Alexander Moroz