      PROGRAM MIECOLOR 
c computes colour of colloidal metal solution
	implicit none
C Parameters:
      INTEGER MXNANG
      PARAMETER(MXNANG=2)
C Variables:
      INTEGER i,J,ju,lambi,NANG,nstuetz
      REAL DANG,GSCA,PI,
     &     QBACK,QEXT,QSCA,RAD,REFMED,
     &     WAVEL,X,brere,breim,ximax,
     &     wlstue(25),qextar(95),
     &     rbar(95),gbar(95),bbar(95)
      COMPLEX REFREL,S1(2*MXNANG-1),S2(2*MXNANG-1),
     &     refind(25)
	character*12 expdatsrc,imagefile
      common /colmat/ rbar,gbar,bbar
C***********************************************************************
      OPEN(UNIT=5,FILE='miecolin.txt',STATUS='UNKNOWN')
	read(5,*) expdatsrc,imagefile  ! data set names, see next two lines
	open(unit=4,file=expdatsrc,status='old')
      OPEN(UNIT=8,FILE=imagefile,STATUS='UNKNOWN')
      READ(5,*)REFMED,ximax
      READ(5,*)RAD
      PI=4.E0*ATAN(1.E0)
      NANG=2
C**********
C NANG=number of angles between 0 and 90 degrees (incl. 0 and 90)
C Scattering matrix elements are calculated for 2*NANG-1 angles
C including 0, 90, and 180 degrees.
C**********
      IF(NANG.GT.1)DANG=0.5E0*PI/FLOAT(NANG-1)

	read(4,*)nstuetz
	do i=1,nstuetz
	read(4,*)wlstue(i),brere,breim
	refind(i)=cmplx(brere,breim)
	enddo
	do lambi=1,95
	wavel=0.355+lambi*0.005
	do j=1,nstuetz
	if(wlstue(j).le.wavel) ju=j
	enddo
	refrel=(refind(ju)*(wlstue(ju+1)-wavel)
     &        +refind(ju+1)*(wavel-wlstue(ju)))/
     &      (wlstue(ju+1)-wlstue(ju))
	refrel=refrel/refmed
      x=2.e0*pi*rad*refmed/wavel

      CALL BHMIE(X,REFREL,NANG,S1,S2,QEXT,QSCA,QBACK,GSCA)

	qextar(lambi)=qext
	enddo

c    create PostScript image:
      write(8,*)'%!PS-Adobe-3.0 EPSF-3.0'
      write(8,*)'%%BoundingBox: 0 0 2048 384'
      write(8,*)'gsave'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'/Hintergrundfarbe{0.133 dup dup}def'
      write(8,*)'/Schriftfarbe {0.7 dup dup} def'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'% Gamma-Korrektur'
      write(8,*)'/Gamma 2.2 def'
      write(8,*)'/YnachL {1 Gamma div exp} def'
      write(8,*)'              % Screen-luminosity-L as a function'
      write(8,*)'              % of the colorimetric brightness Y'
      write(8,*)'/LnachY {Gamma exp} def'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'% Tables of sRGB colour matching functions,'
      write(8,*)'% 5 nm intervals, starting at 360 nm'
      write(8,*)'/rbar 95 array def'
      write(8,8010) (rbar(i),i=1,95)
      write(8,*)'rbar astore pop'
      write(8,*)'/gbar 95 array def'
      write(8,8010) (gbar(i),i=1,95)
      write(8,*)'gbar astore pop'
      write(8,*)'/bbar 95 array def'
      write(8,8010) (bbar(i),i=1,95)
      write(8,*)'bbar astore pop'
 8010 format(5e14.6)
      write(8,*)'/Qext 95 array def'
      write(8,8010) (qextar(i),i=1,95)
      write(8,*)'Qext astore pop'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'% Constants for black-body radiation,'
      write(8,*)'% (for function Slamb in lambInt)'
      write(8,*)'/alf 14.403e6 def'
      write(8,*)'/lamb0 560 def'
      write(8,*)'/Tabs 6504 def'
      write(8,*)'/e 2.71828183 def'
      write(8,*)'/Norm e alf Tabs lamb0 mul div exp 1 sub def'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'/Transmiss{Qext i get'
      write(8,*)'xi mul neg e exch exp'
      write(8,*)'}def'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'/lambInt{'
      write(8,*)'/Rw 0 def /Gw 0 def /Bw 0 def'
      write(8,*)'0 1 94{/i exch def'
      write(8,*)'/lamb 360 5 i mul add def'
      write(8,*)'/Slamb lamb0 lamb div 5 exp Norm mul'
      write(8,*)'e alf Tabs lamb mul div exp 1 sub div def'
      write(8,*)'Transmiss dup dup'
      write(8,*)'rbar i get mul Slamb mul Rw add /Rw exch def'
      write(8,*)'gbar i get mul Slamb mul Gw add /Gw exch def'
      write(8,*)'bbar i get mul Slamb mul Bw add /Bw exch def'
      write(8,*)'}for'
      write(8,*)'}def'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'/lambInt0{'
      write(8,*)'/Rw 0 def /Gw 0 def /Bw 0 def'
      write(8,*)'0 1 94{/i exch def'
      write(8,*)'/lamb 360 5 i mul add def'
      write(8,*)'/Slamb lamb0 lamb div 5 exp'       
      write(8,*)'Norm mul e alf Tabs lamb mul div exp 1 sub div def'
      write(8,*)'1 dup dup'
      write(8,*)'rbar i get mul Slamb mul Rw add /Rw exch def'
      write(8,*)'gbar i get mul Slamb mul Gw add /Gw exch def'
      write(8,*)'bbar i get mul Slamb mul Bw add /Bw exch def'
      write(8,*)'}for'
      write(8,*)'}def'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'Hintergrundfarbe setrgbcolor'
      write(8,*)'0 0 2048 384 rectfill'
      write(8,*)'10 110 translate'
      write(8,*)'/dasMaximum 0 def'
      write(8,*)'lambInt0'
      write(8,*)'Rw dasMaximum gt{/dasMaximum Rw def}if'
      write(8,*)'Gw dasMaximum gt{/dasMaximum Gw def}if'
      write(8,*)'Bw dasMaximum gt{/dasMaximum Bw def}if'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'/ximax ',ximax,' def'
      write(8,*)'/Skalf 2050 ximax div def'
      write(8,*)'0 0.005 ximax {/xi exch def '
      write(8,*)'lambInt'
      write(8,*)'/Rw Rw dasMaximum div def'
      write(8,*)'/Gw Gw dasMaximum div def'
      write(8,*)'/Bw Bw dasMaximum div def'
      write(8,*)'Rw 0 lt {/Rw 0 def} if'
      write(8,*)'Gw 0 lt {/Gw 0 def} if'
      write(8,*)'Bw 0 lt {/Bw 0 def} if'
      write(8,*)'Rw YnachL Gw YnachL Bw YnachL setrgbcolor'
      write(8,*)'/xkoord xi  Skalf mul 1 sub def'
      write(8,*)'xkoord 0 4 270 rectfill'
      write(8,*)'}for'
      write(8,*)'%----------------------------------------------------'
      write(8,*)'/Helvetica findfont 40 scalefont setfont'
      write(8,*)'/strxi 4 string def'
      write(8,*)'0 -50 translate'
      write(8,*)'Schriftfarbe setrgbcolor'
      write(8,*)'0 0.5 17 {/xi exch def'
      write(8,*)'/xkoord xi Skalf mul 2 sub def'
      write(8,*)'xkoord 0 4 45 rectfill'
      write(8,*)'}for'
      write(8,*)'0 0.1 17 {/xi exch def /xkoord xi'
      write(8,*)'Skalf mul 2 sub def xkoord 15 4 30 rectfill}for'
      write(8,*)'0 1 8{/xi exch def'
      write(8,*)'/xkoord xi Skalf mul 2 sub def'
      write(8,*)'xkoord -20 4 30 rectfill'
      write(8,*)'}for'
      write(8,*)'0 1 17  {/xi exch def'
      write(8,*)'/xkoord xi Skalf mul 1 sub def'
      write(8,*)'xkoord 43 sub -55 moveto'
      write(8,*)'xi 10 lt{23 0 rmoveto}if'
      write(8,*)'xi 999 lt{10 0 rmoveto}if'
      write(8,*)'xi strxi cvs show'
      write(8,*)'}for'
      write(8,*)'showpage'
      write(8,*)'grestore'
      write(8,*)'%%EOF'

      STOP
      END
c ---------------------------------------------------------------------
      SUBROUTINE BHMIE(X,REFREL,NANG,S1,S2,QEXT,QSCA,QBACK,GSCA)
	implicit none
C Declare parameters:
C Note: important that MXNANG be consistent with dimension of S1 and S2
C       in calling routine!
      INTEGER MXNANG,NMXX
      PARAMETER(MXNANG=2,NMXX=15000) ! change these parameters if necessary
c      PARAMETER(MXNANG=1000,NMXX=150000)
C Arguments:
      INTEGER NANG
      REAL GSCA,QBACK,QEXT,QSCA,X
      COMPLEX REFREL
      COMPLEX S1(2*MXNANG-1),S2(2*MXNANG-1)
C Local variables:
      INTEGER J,JJ,N,NSTOP,NMX,NN
      DOUBLE PRECISION CHI,CHI0,CHI1,DANG,DX,EN,FN,P,PII,PSI,PSI0,PSI1,
     &                 THETA,XSTOP,YMOD
      DOUBLE PRECISION AMU(MXNANG),PI(MXNANG),PI0(MXNANG),PI1(MXNANG),
     &                 TAU(MXNANG)
      DOUBLE COMPLEX AN,AN1,BN,BN1,DREFRL,XI,XI1,Y
      DOUBLE COMPLEX D(NMXX)
C***********************************************************************
C Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine
C    to calculate scattering and absorption by a homogenous isotropic
C    sphere.
C Given:
C    X = 2*pi*a/lambda
C    REFREL = (complex refr. index of sphere)/(real index of medium)
C    NANG = number of angles between 0 and 90 degrees
C           (will calculate 2*NANG-1 directions from 0 to 180 deg.)
C           if called with NANG<2, will set NANG=2 and will compute
C           scattering for theta=0,90,180.
C Returns:
C    S1(1 - 2*NANG-1) = -i*f_22 (incid. E perp. to scatt. plane,
C                                scatt. E perp. to scatt. plane)
C    S2(1 - 2*NANG-1) = -i*f_11 (incid. E parr. to scatt. plane,
C                                scatt. E parr. to scatt. plane)
C    QEXT = C_ext/pi*a**2 = efficiency factor for extinction
C    QSCA = C_sca/pi*a**2 = efficiency factor for scattering
C    QBACK = (dC_sca/domega)/pi*a**2
C          = backscattering efficiency [NB: this is (1/4*pi) smaller
C            than the "radar backscattering efficiency"; see Bohren &
C            Huffman 1983 pp. 120-123]
C    GSCA = <cos(theta)> for scattering
C
C Original program taken from Bohren and Huffman (1983), Appendix A
C Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26
C in order to compute <cos(theta)>
C 91/05/07 (BTD): Modified to allow NANG=1
C 91/08/15 (BTD): Corrected error (failure to initialize P)
C 91/08/15 (BTD): Modified to enhance vectorizability.
C 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1
C 91/08/15 (BTD): Changed definition of QBACK.
C 92/01/08 (BTD): Converted to full double precision and double complex
C                 eliminated 2 unneed lines of code
C                 eliminated redundant variables (e.g. APSI,APSI0)
C                 renamed RN -> EN = double precision N
C                 Note that DOUBLE COMPLEX and DCMPLX are not part
C                 of f77 standard, so this version may not be fully
C                 portable.  In event that portable version is
C                 needed, use src/bhmie_f77.f
C 93/06/01 (BTD): Changed AMAX1 to generic function MAX
C***********************************************************************
C*** Safety checks
      IF(NANG.GT.MXNANG)STOP'***Error: NANG > MXNANG in bhmie'
      IF(NANG.LT.2)NANG=2
C*** Obtain pi:
      PII=4.*ATAN(1.D0)
      DX=X
      DREFRL=REFREL
      Y=X*DREFRL
      YMOD=ABS(Y)
C
C*** Series expansion terminated after NSTOP terms
C    Logarithmic derivatives calculated from NMX on down
      XSTOP=X+4.*X**0.3333+2.
      NMX=MAX(XSTOP,YMOD)+15
C BTD experiment 91/1/15: add one more term to series and compare results
C      NMX=AMAX1(XSTOP,YMOD)+16
C test: compute 7001 wavelengths between .0001 and 1000 micron
C for a=1.0micron SiC grain.  When NMX increased by 1, only a single
C computed number changed (out of 4*7001) and it only changed by 1/8387
C conclusion: we are indeed retaining enough terms in series!
      NSTOP=XSTOP
C
      IF(NMX.GT.NMXX)THEN
          WRITE(0,*)'Error: NMX > NMXX=',NMXX,' for |m|x=',YMOD
          STOP
      ENDIF
C*** Require NANG.GE.1 in order to calculate scattering intensities
      DANG=0.
      IF(NANG.GT.1)DANG=.5*PII/DBLE(NANG-1)
      DO 1000 J=1,NANG
          THETA=DBLE(J-1)*DANG
          AMU(J)=COS(THETA)
 1000 CONTINUE
      DO 1100 J=1,NANG
          PI0(J)=0.
          PI1(J)=1.
 1100 CONTINUE
      NN=2*NANG-1
      DO 1200 J=1,NN
          S1(J)=(0.,0.)
          S2(J)=(0.,0.)
 1200 CONTINUE
C
C*** Logarithmic derivative D(J) calculated by downward recurrence
C    beginning with initial value (0.,0.) at J=NMX
C
      D(NMX)=(0.,0.)
      NN=NMX-1
      DO 2000 N=1,NN
          EN=NMX-N+1
          D(NMX-N)=(EN/Y)-(1./(D(NMX-N+1)+EN/Y))
 2000 CONTINUE
C
C*** Riccati-Bessel functions with real argument X
C    calculated by upward recurrence
C
      PSI0=COS(DX)
      PSI1=SIN(DX)
      CHI0=-SIN(DX)
      CHI1=COS(DX)
      XI1=DCMPLX(PSI1,-CHI1)
      QSCA=0.E0
      GSCA=0.E0
      P=-1.
      DO 3000 N=1,NSTOP
          EN=N
          FN=(2.E0*EN+1.)/(EN*(EN+1.))
C for given N, PSI  = psi_n        CHI  = chi_n
C              PSI1 = psi_{n-1}    CHI1 = chi_{n-1}
C              PSI0 = psi_{n-2}    CHI0 = chi_{n-2}
C Calculate psi_n and chi_n
          PSI=(2.E0*EN-1.)*PSI1/DX-PSI0
          CHI=(2.E0*EN-1.)*CHI1/DX-CHI0
          XI=DCMPLX(PSI,-CHI)
C
C*** Store previous values of AN and BN for use
C    in computation of g=<cos(theta)>
          IF(N.GT.1)THEN
              AN1=AN
              BN1=BN
          ENDIF
C
C*** Compute AN and BN:
          AN=(D(N)/DREFRL+EN/DX)*PSI-PSI1
          AN=AN/((D(N)/DREFRL+EN/DX)*XI-XI1)
          BN=(DREFRL*D(N)+EN/DX)*PSI-PSI1
          BN=BN/((DREFRL*D(N)+EN/DX)*XI-XI1)
C
C*** Augment sums for Qsca and g=<cos(theta)>
          QSCA=QSCA+(2.*EN+1.)*(ABS(AN)**2+ABS(BN)**2)
          GSCA=GSCA+((2.*EN+1.)/(EN*(EN+1.)))*
     &         (REAL(AN)*REAL(BN)+IMAG(AN)*IMAG(BN))
          IF(N.GT.1)THEN
              GSCA=GSCA+((EN-1.)*(EN+1.)/EN)*
     &        (REAL(AN1)*REAL(AN)+IMAG(AN1)*IMAG(AN)+
     &         REAL(BN1)*REAL(BN)+IMAG(BN1)*IMAG(BN))
          ENDIF
C
C*** Now calculate scattering intensity pattern
C    First do angles from 0 to 90
          DO 2500 J=1,NANG
              JJ=2*NANG-J
              PI(J)=PI1(J)
              TAU(J)=EN*AMU(J)*PI(J)-(EN+1.)*PI0(J)
              S1(J)=S1(J)+FN*(AN*PI(J)+BN*TAU(J))
              S2(J)=S2(J)+FN*(AN*TAU(J)+BN*PI(J))
 2500     CONTINUE
C
C*** Now do angles greater than 90 using PI and TAU from
C    angles less than 90.
C    P=1 for N=1,3,...; P=-1 for N=2,4,...
          P=-P
          DO 2600 J=1,NANG-1
              JJ=2*NANG-J
              S1(JJ)=S1(JJ)+FN*P*(AN*PI(J)-BN*TAU(J))
              S2(JJ)=S2(JJ)+FN*P*(BN*PI(J)-AN*TAU(J))
 2600     CONTINUE
          PSI0=PSI1
          PSI1=PSI
          CHI0=CHI1
          CHI1=CHI
          XI1=DCMPLX(PSI1,-CHI1)
C
C*** Compute pi_n for next value of n
C    For each angle J, compute pi_n+1
C    from PI = pi_n , PI0 = pi_n-1
          DO 2800 J=1,NANG
              PI1(J)=((2.*EN+1.)*AMU(J)*PI(J)-(EN+1.)*PI0(J))/EN
              PI0(J)=PI(J)
 2800     CONTINUE
 3000 CONTINUE
C
C*** Have summed sufficient terms.
C    Now compute QSCA,QEXT,QBACK,and GSCA
      GSCA=2.*GSCA/QSCA
      QSCA=(2./(DX*DX))*QSCA
      QEXT=(4./(DX*DX))*REAL(S1(1))
      QBACK=(ABS(S1(2*NANG-1))/DX)**2/PII
      RETURN
      END
c ---------------------------------------------------------------------
      block data
c       Table of colour matching functions, 5 nm intervals, starting at 360 nm
c       Source:  http://cvrl.ioo.ucl.ac.uk/database/data/cmfs/ciexyz31.txt
c       Sums normalized to 1.0, transformed to sRGB primaries
      real rbar(95),gbar(95),bbar(95)
      common /colmat/ rbar,gbar,bbar
        data (rbar(i),i=1,95)/ 
     *   0.528324E-05,  0.937118E-05,  0.166472E-04,  0.295850E-04,
     *   0.542359E-04,  0.884560E-04,  0.167246E-03,  0.300086E-03,
     *   0.559337E-03,  0.900829E-03,  0.167458E-02,  0.295686E-02,
     *   0.503522E-02,  0.781333E-02,  0.990537E-02,  0.107563E-01,
     *   0.104178E-01,  0.906551E-02,  0.692056E-02,  0.419909E-02,
     *   0.850218E-03, -0.288128E-02, -0.695657E-02, -0.108591E-01,
     *  -0.144672E-01, -0.177740E-01, -0.209726E-01, -0.246289E-01,
     *  -0.288542E-01, -0.339083E-01, -0.384873E-01, -0.419703E-01,
     *  -0.433327E-01, -0.418016E-01, -0.379201E-01, -0.322949E-01,
     *  -0.250810E-01, -0.162980E-01, -0.605381E-02,  0.557875E-02,
     *   0.184923E-01,  0.324298E-01,  0.470596E-01,  0.619003E-01,
     *   0.763697E-01,  0.896935E-01,  0.101207E+00,  0.110292E+00,
     *   0.115731E+00,  0.117838E+00,  0.115912E+00,  0.110621E+00,
     *   0.102218E+00,  0.909029E-01,  0.783962E-01,  0.666025E-01,
     *   0.553634E-01,  0.447967E-01,  0.353139E-01,  0.273101E-01,
     *   0.206297E-01,  0.151810E-01,  0.109579E-01,  0.798018E-02,
     *   0.587280E-02,  0.413394E-02,  0.285334E-02,  0.199148E-02,
     *   0.142828E-02,  0.101985E-02,  0.728066E-03,  0.516714E-03,
     *   0.364556E-03,  0.257661E-03,  0.181059E-03,  0.125731E-03,
     *   0.867691E-04,  0.598538E-04,  0.417828E-04,  0.295266E-04,
     *   0.208914E-04,  0.147632E-04,  0.104457E-04,  0.738162E-05,
     *   0.521936E-05,  0.369081E-05,  0.259948E-05,  0.183071E-05,
     *   0.128931E-05,  0.908009E-06,  0.639486E-06,  0.450348E-06,
     *   0.317176E-06,  0.223374E-06,  0.157316E-06/
        data (gbar(i),i=1,95)/ 
     *  -0.436774E-05, -0.780112E-05, -0.139415E-04, -0.249153E-04,
     *  -0.460647E-04, -0.752557E-04, -0.142872E-03, -0.257421E-03,
     *  -0.482172E-03, -0.781059E-03, -0.146340E-02, -0.260668E-02,
     *  -0.448681E-02, -0.707729E-02, -0.916062E-02, -0.102616E-01,
     *  -0.103766E-01, -0.970079E-02, -0.846387E-02, -0.684733E-02,
     *  -0.467496E-02, -0.192942E-02,  0.162930E-02,  0.546372E-02,
     *   0.944372E-02,  0.134274E-01,  0.177080E-01,  0.227141E-01,
     *   0.286520E-01,  0.360473E-01,  0.440276E-01,  0.522716E-01,
     *   0.595912E-01,  0.647512E-01,  0.682253E-01,  0.701077E-01,
     *   0.705938E-01,  0.697465E-01,  0.676795E-01,  0.645534E-01,
     *   0.603728E-01,  0.551277E-01,  0.489973E-01,  0.421397E-01,
     *   0.348094E-01,  0.272719E-01,  0.199048E-01,  0.130766E-01,
     *   0.721904E-02,  0.233729E-02, -0.131259E-02, -0.382557E-02,
     *  -0.530242E-02, -0.589611E-02, -0.586896E-02, -0.552474E-02,
     *  -0.494899E-02, -0.422962E-02, -0.346307E-02, -0.275429E-02,
     *  -0.212292E-02, -0.158267E-02, -0.115423E-02, -0.847488E-03,
     *  -0.628550E-03, -0.445523E-03, -0.308668E-03, -0.215907E-03,
     *  -0.155012E-03, -0.110686E-03, -0.790180E-04, -0.560796E-04,
     *  -0.395658E-04, -0.279643E-04, -0.196505E-04, -0.136458E-04,
     *  -0.941725E-05, -0.649598E-05, -0.453474E-05, -0.320458E-05,
     *  -0.226737E-05, -0.160227E-05, -0.113370E-05, -0.801135E-06,
     *  -0.566462E-06, -0.400568E-06, -0.282126E-06, -0.198690E-06,
     *  -0.139931E-06, -0.985471E-07, -0.694044E-07, -0.488768E-07,
     *  -0.344235E-07, -0.242431E-07, -0.170736E-07/
        data (bbar(i),i=1,95)/ 
     *   0.302874E-04,  0.542673E-04,  0.972396E-04,  0.174190E-03,
     *   0.322301E-03,  0.527168E-03,  0.100187E-02,  0.180932E-02,
     *   0.339034E-02,  0.550637E-02,  0.103628E-01,  0.185513E-01,
     *   0.322525E-01,  0.518961E-01,  0.691805E-01,  0.809898E-01,
     *   0.871223E-01,  0.888151E-01,  0.881869E-01,  0.866600E-01,
     *   0.827671E-01,  0.755501E-01,  0.633452E-01,  0.508421E-01,
     *   0.391416E-01,  0.290201E-01,  0.211111E-01,  0.150477E-01,
     *   0.103849E-01,  0.661958E-02,  0.304696E-02, -0.206790E-03,
     *  -0.274501E-02, -0.445793E-02, -0.571602E-02, -0.667338E-02,
     *  -0.735090E-02, -0.776296E-02, -0.794085E-02, -0.793279E-02,
     *  -0.776180E-02, -0.744355E-02, -0.700364E-02, -0.645951E-02,
     *  -0.584110E-02, -0.517832E-02, -0.450256E-02, -0.383521E-02,
     *  -0.322046E-02, -0.266043E-02, -0.217594E-02, -0.175785E-02,
     *  -0.140403E-02, -0.110382E-02, -0.855329E-03, -0.659655E-03,
     *  -0.503830E-03, -0.379709E-03, -0.283550E-03, -0.209732E-03,
     *  -0.153110E-03, -0.110104E-03, -0.779897E-04, -0.559292E-04,
     *  -0.405493E-04, -0.281566E-04, -0.192890E-04, -0.134034E-04,
     *  -0.959207E-05, -0.684905E-05, -0.488954E-05, -0.347015E-05,
     *  -0.244827E-05, -0.173040E-05, -0.121596E-05, -0.844387E-06,
     *  -0.582714E-06, -0.401969E-06, -0.280605E-06, -0.198293E-06,
     *  -0.140303E-06, -0.991476E-07, -0.701503E-07, -0.495738E-07,
     *  -0.350524E-07, -0.247869E-07, -0.174576E-07, -0.122947E-07,
     *  -0.865873E-08, -0.609806E-08, -0.429465E-08, -0.302446E-08,
     *  -0.213011E-08, -0.150013E-08, -0.105651E-08/
      end
