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