c     Mie-colours: light scattering by a spherical water drop, 
c     incident light: blackbody radiation. 
c	Output is an encapsulated PostScript image.
c
c
      implicit real*8 (a-h,o-z)
      parameter(NMXX=1500,lmxx=200,nmxradii=31)
      real*8 Gwicht(50),sl(lmxx),Ru(nmxradii),weight(nmxradii)
      real*4 winklf(3,NMXX)
	real*4 s1qps2q(NMXX),sqave(NMXX,lmxx)
	complex*16 H2On
	logical Test,lognorm
	integer plotscale
	character*12 imagename
      common XuR,YuR,ZuR,XuG,YuG,ZuG,XuB,YuB,ZuB,
     1        RuX,GuX,BuX,RuY,GuY,BuY,RuZ,GuZ,BuZ,
     2        Ru,weight,pi,grad,
     3        Gwicht,ntotlamb

      open(unit=5,file='rainbdat.txt',status='unknown')
	read(5,*) imagename
      read(5,*) Tabs,nstart,nstop,Nprograd,ntotlamb,plotscale
c Tabs = colour temperature, nstart, nstop: angles in degrees from antisolar point
c Nprograd: number of angles per degree, ntotlamb: number of wavelengths
c plotscale: number of pixels per degree
      read(5,*) grundR,grundG,grundB,lognorm,Test
c grundR,grundG,grundB: background colour
c lognorm: logical variable. If lognorm, then Gaussian distribution of the 
c   logarithms of droplet radius is assumed, else Gaussian distribution of radii.
c Test: logical variable. If Test, only the radii and weights are written to
c   the screen.
	read(5,*) numRad,Rmostp,Faktor,sigma
c numRad = number of bows to be superposed, Rmostp = most probable radius
c Faktor = (most probable radius)/(second frequent radius)
c sigma = approximate standard deviation (if numRad is large enough)
	if(numrad.gt.nmxradii)then
	write(*,*) 'numrad too large or nmxradii too small'
	stop
	endif
	Summand=(Faktor-1.0d0)*Rmostp
	sigmalog=sigma/Rmostp
	write(*,*)'   radii and relative weights:'
	if (lognorm) then
	do i=1,numRad
	Ru(i)=Rmostp*Faktor**(real(i)-real(numRad+1)/2.0)
	weight(i)=exp(-(log(Ru(i)/Rmostp))**2/(2.0*sigmalog**2))
	write(*,*)Ru(i),weight(i)
	enddo
	else
	do i=1,numRad
	Ru(i)=Rmostp+Summand*(real(i)-real(numRad+1)/2.0)
	weight(i)=exp(-(Ru(i)-Rmostp)**2/(2.0*sigma**2))
	write(*,*)Ru(i),weight(i)
	enddo
	endif
	Rave=0.0
	Rsqave=0.0
	weignor=0.0
	do i=1,numRad
	Rave=Rave+Ru(i)*weight(i)
	Rsqave=Rsqave+Ru(i)**2*weight(i)
	weignor=weignor+weight(i)
	enddo
	Rave=Rave/weignor
	Rsqave=Rsqave/weignor
	deltaR=sqrt(Rsqave-Rave**2)
	write(*,*)'Rave=',Rave,'   deltaR=',deltaR	
	if (Test) then
	write(*,*) 'test of size distribution only'
	stop
	endif

      open(unit=8,file=imagename,status='unknown')

      SonnRad=0.267d0
	maddang=int(SonnRad*Nprograd)
	ndifferenz=nstop-nstart
	nang=ndifferenz*Nprograd+2*maddang
c
       pi=4*atan(1.0d0)
       grad=pi/180.0d0

c	normalized blackbody radiation:
      alf=14.388d3
      alamo=0.560d0
        do i=1,ntotlamb
        wlambda=0.360d0+0.47d0*real(i-1)/real(ntotlamb-1)
        sl(i)=(alamo/wlambda)**5*(exp(alf/(Tabs*alamo))-1.0d0)
     +      /(exp(alf/(Tabs*wlambda))-1.0d0)
	enddo

c	sRGB primaries
	xr=0.64
	yr=0.33
	xg=0.30
	yg=0.60
	xb=0.15
	yb=0.06
c       white D65
       xw=0.3127
       yw=0.3291

c
        Gamma=2.2
       zr=1.d0-(xr+yr)
       zg=1.d0-(xg+yg)
       zb=1.d0-(xb+yb)

       zw=1.d0-(xw+yw)

c       calibrate the primaries to white
c       (As FORTRAN does not discriminate between upper- and lowercase,
c       instead of X, Xu is written, u means uppercase)
        XuW=xw/yw
        YuW=1.0d0
        ZuW=zw/yw
        call LG(xr,xg,xb,yr,yg,yb,zr,zg,zb,XuW,YuW,ZuW,u1,u2,u3)
        XuR=xr*u1
        YuR=yr*u1
        ZuR=zr*u1
        XuG=xg*u2
        YuG=yg*u2
        ZuG=zg*u2
        XuB=xb*u3
        YuB=yb*u3
        ZuB=zb*u3

c       Determination of RGB-coordinates of the virtual primaries XYZ.
c       Matrix inversion done in three steps.

c       superpose X from R, G, B
        call LG(XuR,XuG,XuB,YuR,YuG,YuB,ZuR,ZuG,ZuB,1.d0,0.d0,0.d0,
     1     RuX,GuX,BuX)
c       superpose Y from R, G, B
        call LG(XuR,XuG,XuB,YuR,YuG,YuB,ZuR,ZuG,ZuB,0.d0,1.d0,0.d0,
     1     RuY,GuY,BuY)
c       superpose Z from R, G, B
        call LG(XuR,XuG,XuB,YuR,YuG,YuB,ZuR,ZuG,ZuB,0.d0,0.d0,1.d0,
     1     RuZ,GuZ,BuZ)

c       transformation matrix ready.
c       Transformation from XYZ to RGB is
c       RuColour = RuX*XuColour + RuY*YuColour + RuZ*ZuColour
c       and accordingly for G and B

c------------------------------------------------------------------------------

	  hinterR=grundR**Gamma
	  hinterG=grundG**Gamma
	  hinterB=grundB**Gamma
        nhor=ndifferenz*Nprograd
	if(nhor.gt.NMXX-2*maddang)then
	write(*,*) 'NHOR is too large. Reduce number of angles'
	stop
	endif

c     Prepare averaging over the sun's disc
      Schritt=1.d0/Nprograd      ! step (degrees) 
      do i=1, 2*maddang+1
      Gwicht(i)=sqrt(max(0.d0,SonnRad**2-((i-maddang-1)*Schritt)**2))
      enddo
      dnorm=0.0d0
      do i=1, 2*maddang+1
      dnorm=dnorm+Gwicht(i)
      enddo
      dnorm=1.d0/dnorm
      do i=1, 2*maddang+1
      Gwicht(i)=Gwicht(i)*dnorm
      enddo

	astart=(dble(nstart)-maddang*Schritt)*grad
	astop=(dble(nstop)+maddang*Schritt)*grad

c  	get angular distribution for all wavelengths
	almstep=0.470d0/(ntotlamb-1)
	do ilamb=1,ntotlamb
	wlambda=0.360d0+(ilamb-1)*almstep
      H2On=dcmplx(1.3175+0.007382/(wlambda-0.11278))
	do j=1, nang-2*maddang
	sqave(j,ilamb)=0.0
	enddo
	do m=1,numRad
      wkRu=2*pi*Ru(m)/wlambda
	call BHMIEMOD(wkRu,H2On,NANG,astart,astop,wlambda,s1qps2q)
	do j=1, nang-2*maddang
	do jj=1,2*maddang+1
      sqave(j,ilamb)=sqave(j,ilamb)
     &     +s1qps2q(j+jj-1)*real(Gwicht(jj))*weight(m)
	enddo !jj
	enddo !j
	enddo !m
	enddo !ilamb

	

c       compute colour coordinates for all angles
c	  and store them in array winklf
        do iw=1,nang-2*maddang
        call lambInt(XF,YF,ZF,iw,sl,sqave)
c       XYZ to RGB
         RuF=RuX*XF+RuY*YF+RuZ*ZF
         GuF=GuX*XF+GuY*YF+GuZ*ZF
         BuF=BuX*XF+BuY*YF+BuZ*ZF
           winklf(1,iw)=real(RuF)
           winklf(2,iw)=real(GuF)
           winklf(3,iw)=real(BuF)
        enddo   ! iw    


        dasRmaximum=0.d0
        dasGmaximum=0.d0
        dasBmaximum=0.d0
        do j=1,nhor
        if (winklf(1,j) .gt. dasRmaximum) dasRmaximum=winklf(1,j) 
        if (winklf(2,j) .gt. dasGmaximum) dasGmaximum=winklf(2,j) 
        if (winklf(3,j) .gt. dasBmaximum) dasBmaximum=winklf(3,j) 
        enddo

c       Create PostScript-image
      write(8,*) '%!PS-Adobe-3.0 EPSF-3.0 '
      write(8,*) '%%BoundingBox: 0 0 4100 410'
      write(8,*) '% Mie scattering intensity computed by MieColor.for'
      write(8,1) Rmostp,sigma
    1 format('% water drops, radius Rmax=',f9.3,',  sigma=',
     *     f9.3,' (microns)') 
      write(8,*) 'gsave '
      write(8,*)'/nhor ',nhor,' def /Nprograd ',Nprograd,' def'
      write(8,*) '/plotscale ',plotscale,' def'
      write(8,*)'%----------------------------------------------------'
      write(8,*) '/Gamma ',Gamma,' def          % Gamma of the screen '
      write(8,*) '                        % may be adjusted!        '                      
      write(8,*) '% Gamma-Korrektur '
      write(8,*) '/YnachL {1 Gamma div exp} def '
      write(8,*) '      % Brightness on the screen L is a function '
      write(8,*) '      % of the colorimetric brightness Y'
      write(8,*) '/LnachY {Gamma exp} def'
      write(8,*)'%----------------------------------------------------'
      write(8,*) '% Background colour:'
      write(8,3) grundR,grundG,grundB
    3 format ('/hinterR ',f7.5,' def /hinterG ',
     &    f7.5,' def /hinterB ', f7.5,' def')
      write(8,*) '0 0 4100 410 rectclip'
      write(8,*)
     +   '0.1333 dup dup setrgbcolor 0 0 4100 410 rectfill'
	write(8,*) '/hinterR hinterR LnachY def'
	write(8,*) '/hinterG hinterG LnachY def'
	write(8,*) '/hinterB hinterB LnachY def'
    4 format('/dasRmaximum ',e12.5,' def /dasGmaximum ',
     &    e12.5,' def /dasBmaximum ',e12.5,' def')     
      write(8,4) dasRmaximum,dasGmaximum,dasBmaximum
      write(8,*) '/dasRmaximum dasRmaximum 1 hinterR sub div def'
      write(8,*) '/dasGmaximum dasGmaximum 1 hinterG sub div def'
	write(8,*) '/dasBmaximum dasBmaximum 1 hinterB sub div def'
	write(8,*) '/dasMaximum 0 def dasMaximum dasRmaximum lt'
	write(8,*) '{/dasMaximum dasRmaximum def} if'
      write(8,*) 'dasMaximum dasGmaximum lt'
	write(8,*) '{/dasMaximum dasGmaximum def} if'
	write(8,*) 'dasMaximum dasBmaximum lt'
	write(8,*) '{/dasMaximum dasBmaximum def} if'
      write(8,*)'%----------------------------------------------------'
      write(8,*) '50 120 translate'
      write(8,*) '% 0 0 0 setrgbcolor -62 -120 4096 120 rectfill'
      write(8,*) '/pi ', pi,' def'
      write(8,*) -nstart,' plotscale mul 0 translate'
      write(8,*) '/rwert nhor array def'
      write(8,5) (winklf(1,i),i=1,nhor)
    5 format(5e13.5)
      write(8,*) 'rwert astore pop'
      write(8,*) '/gwert nhor array def'
      write(8,5) (winklf(2,i),i=1,nhor)
      write(8,*) 'gwert astore pop'
      write(8,*) '/bwert nhor array def'
      write(8,5) (winklf(3,i),i=1,nhor)
      write(8,*) 'bwert astore pop'
      write(8,*) '/Schritt 1 Nprograd div def'
      write(8,*) '400 Nprograd div setlinewidth'
      write(8,*) nstart,' plotscale mul',
     *     ' currentlinewidth 4 div add 0 moveto'
      write(8,*) '0 1 nhor 1 sub {/i exch def'
      write(8,*) '/Theta Schritt i mul ',nstart,' add def'
      write(8,*)
     +  'rwert i get dasMaximum div hinterR add ',
     +    'dup 0 lt {pop 0} if YnachL'
      write(8,*)
     +   'gwert i get dasMaximum div hinterG add ',
     +    'dup 0 lt {pop 0} if YnachL'
      write(8,*)
     +   'bwert i get dasMaximum div hinterB add ',
     +    'dup 0 lt {pop 0} if YnachL setrgbcolor'
      write(8,*)
     + 'Theta plotscale mul currentlinewidth 4 div add 0 moveto'
      write(8,*) '0 400 rlineto'
      write(8,*) 'stroke '
      write(8,*) '}for'
      write(8,*)
     +   '/Helvetica findfont 50 scalefont setfont'
      write(8,*) '0.7 setgray'
      write(8,*) nstart,' 0.25 ',nstop,' {/Theta exch def'
      write(8,*) 'Theta plotscale mul 1 sub -50 3 47 rectfill'
      write(8,*) '}for'
      write(8,*) nstart,' 1  ',nstop,' {/Theta exch def'
      write(8,*) 'Theta plotscale mul 1 sub -60 3 10 rectfill'
      write(8,*) '}for'
      write(8,*) '/strThet 3 string def'
      write(8,*) nstart,' 1  ',nstop,'  {/Theta exch def'
      write(8,*) 'Theta plotscale mul 25 sub -110 moveto'
      write(8,*) 'Theta strThet cvs show'
      write(8,*) '}for'
      write(8,*) '%showpage '
      write(8,*) 'grestore '
      stop
      end 
c----------------------------------------------------------------------
        subroutine lambInt(Xw,Yw,Zw,iw,sl,sqave)
        implicit real*8 (a-h,o-z)
      parameter(NMXX=1500,lmxx=200,nmxradii=31)
        real*4 xbar(95),ybar(95),zbar(95),sqave(nmxx,lmxx)
	real*8 Gwicht(50),sl(lmxx),Ru(nmxradii),weight(nmxradii)
        common XuR,YuR,ZuR,XuG,YuG,ZuG,XuB,YuB,ZuB,
     1        RuX,GuX,BuX,RuY,GuY,BuY,RuZ,GuZ,BuZ,
     2        Ru,weight,pi,grad,
     3        Gwicht,ntotlamb


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
        data (xbar(i),i=1,95)/
     *   0.607772E-05,  0.108594E-04,  0.194122E-04,  0.346977E-04, 
     *   0.640056E-04,  0.104617E-03,  0.198520E-03,  0.357926E-03, 
     *   0.669532E-03,  0.108501E-02,  0.203573E-02,  0.363213E-02, 
     *   0.628733E-02,  0.100486E-01,  0.132830E-01,  0.153698E-01, 
     *   0.162952E-01,  0.162849E-01,  0.157300E-01,  0.149112E-01, 
     *   0.136059E-01,  0.117484E-01,  0.914044E-02,  0.664853E-02, 
     *   0.447477E-02,  0.271135E-02,  0.149767E-02,  0.687779E-03, 
     *   0.229260E-03,  0.112290E-03,  0.435125E-03,  0.136152E-02, 
     *   0.296026E-02,  0.512793E-02,  0.774336E-02,  0.105623E-01, 
     *   0.135871E-01,  0.168295E-01,  0.202801E-01,  0.239576E-01, 
     *   0.278153E-01,  0.317408E-01,  0.356569E-01,  0.394186E-01, 
     *   0.428716E-01,  0.457864E-01,  0.480182E-01,  0.494405E-01, 
     *   0.496979E-01,  0.489212E-01,  0.469093E-01,  0.439056E-01, 
     *   0.399777E-01,  0.351563E-01,  0.300564E-01,  0.253542E-01, 
     *   0.209562E-01,  0.168810E-01,  0.132643E-01,  0.102325E-01, 
     *   0.771529E-02,  0.567067E-02,  0.408924E-02,  0.297570E-02, 
     *   0.218826E-02,  0.153931E-02,  0.106208E-02,  0.741117E-03, 
     *   0.531469E-03,  0.379491E-03,  0.270917E-03,  0.192272E-03, 
     *   0.135653E-03,  0.958769E-04,  0.673729E-04,  0.467853E-04, 
     *   0.322872E-04,  0.222719E-04,  0.155476E-04,  0.109870E-04, 
     *   0.777380E-05,  0.549348E-05,  0.388690E-05,  0.274674E-05, 
     *   0.194215E-05,  0.137337E-05,  0.967281E-06,  0.681218E-06, 
     *   0.479760E-06,  0.337875E-06,  0.237956E-06,  0.167577E-06, 
     *   0.118023E-06,  0.831187E-07,  0.585380E-07/ 
        data (ybar(i),i=1,95)/
     *   0.183282E-06,  0.325903E-06,  0.579747E-06,  0.103035E-05, 
     *   0.182487E-05,  0.299466E-05,  0.561498E-05,  0.101538E-04, 
     *   0.185294E-04,  0.299466E-04,  0.566177E-04,  0.102005E-03, 
     *   0.187166E-03,  0.341578E-03,  0.542781E-03,  0.787969E-03, 
     *   0.107620E-02,  0.139439E-02,  0.177808E-02,  0.224599E-02, 
     *   0.280749E-02,  0.345789E-02,  0.425709E-02,  0.526872E-02, 
     *   0.650495E-02,  0.792180E-02,  0.973357E-02,  0.121003E-01, 
     *   0.151137E-01,  0.190582E-01,  0.235361E-01,  0.284586E-01, 
     *   0.332220E-01,  0.371150E-01,  0.403343E-01,  0.428072E-01, 
     *   0.446391E-01,  0.458697E-01,  0.465552E-01,  0.467915E-01, 
     *   0.465575E-01,  0.457902E-01,  0.445455E-01,  0.428329E-01, 
     *   0.407086E-01,  0.381959E-01,  0.354212E-01,  0.325154E-01, 
     *   0.295254E-01,  0.265214E-01,  0.235361E-01,  0.206444E-01, 
     *   0.178276E-01,  0.150201E-01,  0.123997E-01,  0.101538E-01, 
     *   0.818851E-02,  0.646658E-02,  0.500669E-02,  0.381819E-02, 
     *   0.285428E-02,  0.208596E-02,  0.149733E-02,  0.108556E-02, 
     *   0.795455E-03,  0.557755E-03,  0.384158E-03,  0.267788E-03, 
     *   0.191939E-03,  0.137052E-03,  0.978410E-04,  0.694386E-04, 
     *   0.489907E-04,  0.346257E-04,  0.243316E-04,  0.168964E-04, 
     *   0.116604E-04,  0.804346E-05,  0.561498E-05,  0.396792E-05, 
     *   0.280749E-05,  0.198396E-05,  0.140374E-05,  0.991980E-06, 
     *   0.701404E-06,  0.495990E-06,  0.349331E-06,  0.246020E-06, 
     *   0.173264E-06,  0.122023E-06,  0.859372E-07,  0.605201E-07, 
     *   0.426238E-07,  0.300181E-07,  0.211409E-07/
        data (zbar(i),i=1,95)/
     *   0.283507E-04,  0.507983E-04,  0.910253E-04,  0.163060E-03, 
     *   0.301703E-03,  0.493482E-03,  0.937851E-03,  0.169374E-02, 
     *   0.317373E-02,  0.515467E-02,  0.970126E-02,  0.173678E-01, 
     *   0.301983E-01,  0.486022E-01,  0.648123E-01,  0.759149E-01, 
     *   0.817198E-01,  0.833822E-01,  0.828915E-01,  0.815813E-01, 
     *   0.780778E-01,  0.714778E-01,  0.602301E-01,  0.487355E-01, 
     *   0.380262E-01,  0.288231E-01,  0.217591E-01,  0.165258E-01, 
     *   0.127230E-01,  0.993046E-02,  0.739990E-02,  0.522484E-02, 
     *   0.366019E-02,  0.267790E-02,  0.197206E-02,  0.139578E-02, 
     *   0.949545E-03,  0.626793E-03,  0.409287E-03,  0.268960E-03, 
     *   0.182425E-03,  0.128633E-03,  0.982288E-04,  0.841961E-04, 
     *   0.771798E-04,  0.654859E-04,  0.514532E-04,  0.467756E-04, 
     *   0.374205E-04,  0.280654E-04,  0.159037E-04,  0.112261E-04, 
     *   0.888737E-05,  0.467756E-05,  0.233878E-05,  0.140327E-05, 
     *   0.935512E-06,  0.467756E-06,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00,  0.000000E+00, 
     *   0.000000E+00,  0.000000E+00,  0.000000E+00/


      Xw=0.0d0   
      Yw=0.0d0   
      Zw=0.0d0 
	diff=0.470d0/(ntotlamb-1)  
	if(ntotlamb.ge.95)then
      do j=1,ntotlamb-1
	almbda=0.360d0+(j-1)*diff
	ku=int((almbda-0.354999999999d0)/0.005d0)
	diffu=(almbda-0.355d0-ku*0.005d0)
	diffo=(0.360d0+ku*0.005d0-almbda)
      hil=dble(sqave(iw,j))*sl(j)
      Xw=Xw+dble(xbar(ku)*diffo+xbar(ku+1)*diffu)*hil
      Yw=Yw+dble(ybar(ku)*diffo+ybar(ku+1)*diffu)*hil
      Zw=Zw+dble(zbar(ku)*diffo+zbar(ku+1)*diffu)*hil
      enddo
	hil=dble(sqave(iw,ntotlamb))*sl(ntotlamb)
	Xw=Xw+xbar(95)*hil*0.005d0
	Yw=Yw+ybar(95)*hil*0.005d0
	Zw=Zw+zbar(95)*hil*0.005d0
	else
	do j=1,94
	almbda=0.355d0+j*0.005d0
	ku=int((almbda-0.359999999999d0+diff)/diff)
	diffu=(almbda-0.360d0+diff-ku*diff)
	diffo=(0.360d0+ku*diff-almbda)
      hil=dble(sqave(iw,ku))*sl(ku)*diffo
     &   +dble(sqave(iw,ku+1))*sl(ku+1)*diffu
      Xw=Xw+xbar(j)*hil
      Yw=Yw+ybar(j)*hil
      Zw=Zw+zbar(j)*hil
	enddo
	hil=sqave(iw,ntotlamb)*sl(ntotlamb)
	Xw=Xw+xbar(95)*hil*diff
	Yw=Yw+ybar(95)*hil*diff
	Zw=Zw+zbar(95)*hil*diff
	endif
      return
      end
c----------------------------------------------------------------------
      subroutine LG(a11,a12,a13,a21,a22,a23,a31,a32,a33,b1,b2,b3,
     1   x,y,z)
      implicit real*8 (a-h,o-z)
      det=a11*a22*a33+a12*a23*a31+a13*a21*a32
     1  -a13*a22*a31-a12*a21*a33-a11*a23*a32
      x =(b1*a22*a33+a12*a23*b3+a13*b2*a32
     1  -a13*a22*b3-a12*b2*a33-b1*a23*a32)/det
      y =(a11*b2*a33+b1*a23*a31+a13*a21*b3
     1  -a13*b2*a31-b1*a21*a33-a11*a23*b3)/det
      z =(a11*a22*b3+a12*b2*a31+b1*a21*a32
     1  -b1*a22*a31-a12*a21*b3-a11*b2*a32)/det
      return
      end
c----------------------------------------------------------------------

      SUBROUTINE BHMIEMOD(X,REFREL,NANG,astart,astop,wlambda,s1qps2q) 
C Declare parameters:
      INTEGER MXNANG,NMXX
C      PARAMETER(MXNANG=1000,NMXX=15000)
      PARAMETER(MXNANG=1500,NMXX=30000)
C Arguments:
      INTEGER NANG
      REAL*8 X,astart,astop,wlambda
	real*4 s1qps2q(mxnang)
      double COMPLEX REFREL
C Local variables:
      double COMPLEX S1(MXNANG),S2(MXNANG)
      INTEGER J,N,NSTOP,NMX,NN
      Real*8 CHI,CHI0,CHI1,DANG,DX,EN,FN,PII,PSI,PSI0,PSI1,
     &                 THETA,XSTOP,YMOD
      Real*8 AMU(MXNANG),PI(MXNANG),PI0(MXNANG),PI1(MXNANG),
     &                 TAU(MXNANG)
      DOUBLE COMPLEX AN,BN,DREFRL,XI,XI1,Y
      DOUBLE COMPLEX D(NMXX)
C***********************************************************************
c    Modification of the 
C Subroutine BHMIE which 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 
C Returns:
C    S1(1 - NANG) = -i*f_22 (incid. E perp. to scatt. plane,
C                                scatt. E perp. to scatt. plane)
C    S2(1 - NANG) = -i*f_11 (incid. E parallel to scatt. plane,
C                                scatt. E parallel to scatt. plane)
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***********************************************************************
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
      NSTOP=XSTOP
C
      IF(NMX.GT.NMXX)THEN
          WRITE(*,*)'Error: NMX > NMXX=',NMXX,' for |m|x=',YMOD
          STOP
      ENDIF
	dang=(astop-astart)/dble(nang-1)
      DO 1000 J=1,NANG
          THETA=astart + DBLE(J-1)*DANG
		theta=pii-theta     !    counting angle from antisolar point
          AMU(J)=COS(THETA)
 1000 CONTINUE
      DO 1100 J=1,NANG
          PI0(J)=0.
          PI1(J)=1.
 1100 CONTINUE
      NN=NANG
      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)
cc      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*** 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*** Now calculate scattering intensity pattern
C    Do angles
          DO 2500 J=1,NANG
              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
          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.
	do 3100 j=1,nang
      s1qps2q(j)=real(ABS(S2(j))*ABS(S2(j))
     &     +ABS(S1(j))*ABS(S1(j)))*wlambda**2
 3100 continue
c   s1qps2q is proportional to scattered intensity. Constant 
c   factors are omitted
      RETURN
      END
