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='glorydat.txt',status='unknown')
	read(5,*) imagename
      read(5,*) Tabs,nstart,nstop,Nprograd,ntotlamb,plotscale		
c Tabs = colour temperature, nstart, nstop: scattering angles in degrees
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.0)*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,'   delta-R=',deltaR	
	if (Test) then
	write(*,*) 'test of size distribution only'
	stop
	endif

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

c	If Test.eq.'true', then colours outside of the sRGB gamut will be
c     marked by false colouring. The background colour (grundR,grundG,grundB)
c	may be adjusted for colour fidelity.

      SonnRad=0.267d0
	maddang=int(SonnRad*Nprograd)
	ndifferenz=nstop-nstart
	nang=ndifferenz*Nprograd+2*maddang
c       Ru ... radius of drop (microns)
c	  Tabs ... colour temperature
c	  Nprograd ... steps per degree
c	nstart, nstop: minimum and maximum angle, measured from the antisolar point

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------------------------------------------------------------------------------
	Gamma=2.2
	  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(RuF,GuF,BuF,iw,sl,sqave)
           winklf(1,iw)=real(RuF)
           winklf(2,iw)=real(GuF)
           winklf(3,iw)=real(BuF)
        enddo   ! iw    


c       Normalize to maximum 1.0,
c       Gamma-correction
        ammaG=1.d0/Gamma
        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 2048 1536'
      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 2.2 def          % Gamma of the screen '
      write(8,*) '                        % may be adjusted!        '                      
      write(8,*) '% Gamma-correction '
      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,*) '1024 768 translate'
      write(8,*) '-1024 -768 2048 1536 rectclip'
      write(8,*)
     +   'hinterR hinterG hinterB setrgbcolor 0 0 4096 768 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,*)'%----------------------------------------------------'
c      write(8,*) '50 120 translate'
      write(8,*) '% 0 0 0 setrgbcolor -62 -120 4096 120 rectfill'
      write(8,*) '/pi ', pi,' def'
      write(8,*) '/tanmul {dup sin exch cos div 180 mul pi div ', 
     &     'plotscale mul} def'
      write(8,*) -nstart,' tanmul 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,' tanmul 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,*)
     +   '0 0 Theta tanmul currentlinewidth 4 div add 0 360 arc'
      write(8,*) 'stroke '
      write(8,*) '}for'
      write(8,*)
     +   '/Helvetica findfont 50 scalefont setfont'
      write(8,*) '0.4 setgray'
      write(8,*) nstart,' 0.25 ',nstop,' {/Theta exch def'
      write(8,*) 'Theta tanmul 1 sub -50 3 47 rectfill'
      write(8,*) '}for'
      write(8,*) nstart,' 1  ',nstop,' {/Theta exch def'
      write(8,*) 'Theta tanmul 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 tanmul 13 sub -110 moveto'
      write(8,*) 'Theta strThet cvs show'
      write(8,*) '}for'
      write(8,*) '%showpage '
      write(8,*) 'grestore '
      stop
      end 
c----------------------------------------------------------------------
        subroutine lambInt(Rw,Gw,Bw,iw,sl,sqave)
        implicit real*8 (a-h,o-z)
      parameter(NMXX=1500,lmxx=200,nmxradii=31)
      real*4 rbar(95),gbar(95),bbar(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, transformed to sRGB primaries
        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/


      Rw=0.0d0   
      Gw=0.0d0   
      Bw=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)
      Rw=Rw+dble(rbar(ku)*diffo
     &     +rbar(ku+1)*diffu)*hil
      Gw=Gw+dble(gbar(ku)*diffo
     &     +gbar(ku+1)*diffu)*hil
      Bw=Bw+dble(bbar(ku)*diffo
     &     +bbar(ku+1)*diffu)*hil
      enddo
	hil=dble(sqave(iw,ntotlamb))*sl(ntotlamb)
	Rw=Rw+rbar(95)*hil*0.005d0
	Gw=Gw+gbar(95)*hil*0.005d0
	Bw=Bw+bbar(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)
      Rw=Rw+rbar(j)*hil
      Gw=Gw+gbar(j)*hil
      Bw=Bw+bbar(j)*hil
	enddo
	hil=sqave(iw,ntotlamb)*sl(ntotlamb)
	Rw=Rw+rbar(95)*hil*diff
	Gw=Gw+gbar(95)*hil*diff
	Bw=Bw+bbar(95)*hil*diff
	endif
      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    wlambda = wavelength
c    astart, astop: first and last angle to compute (radians)
c      (angle measured from the antisolar point)
C Returns:
c    s1qps2q = relative intensity scattered to angle 
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
