C	ARTURO QUIRANTES SIERRA
C	Department of Applied Physics, Faculty of Sciences
C	University of Granada, 18071 Granada (SPAIN)
C	http://www.ugr.es/local/aquiran/codigos.htm
C	aquiran@ugr.es
C
C	Last update: 20 September 2007
C	
C	This is LISA, a computer code to calculate light-scattering properties
C	for randomly-oriented, coated nonspherical scatterers
C	with a non-absorbing (n real) coating
C	It is based on a T-matrix code, combined with Mishchenko's
C	orientation averaging scheme (JOSA A 8, 871-882 (1991)
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, but bugs may remain
C	  and overflow errors, particularly for large sizes and/or extreme
C	  shapes.  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, has you fired,
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	- And a final message.  Don't you think this is the ultimate code.
C	  I've been improving and speeding it up for a long time, and I think
C	  there's still room for improvement.  So, as you grow confidence,
C	  I encourage you to strip it, analize it, and improve it at will.
C	  Improvements, comments, additions: all welcome at aquiran@ugr.es
C
C	  You can find a description of this code, together with benchmark data,
C	  in my JQSRT paper (vol 92, p. 373-381, 2005
C	  available at http://www.ugr.es/~aquiran/ciencia/arti16.pdf
C	  Plus, you'll find plenty of remarks along the code.
C	  Even if you just want to 'plug and play' it's wise to take a look
C	  just to see how you can adapt it to your particular needs.
C
C	And now, let's get down to business.  Enjoy!  Arturo Quirantes.
C
C
	INTEGER SU,SUP,ST,SA,SS,S1,SD,ST1,ST2,SE
	INTEGER SAA,SSS,SUU,SDD,S11,SEE
C	The following parameters set the memory usage
C	SU is the maximum size of the T-matrix (=2*nmax)
C	SA is the number of angles where calculations will be made
C	If you have computer memory problems, I suggest you lower the value
C	of SU.  Remember, the higher SU, the larger the maximum particle size
C	But you have to change SU and/or SA in the entire code, not just here
C	A simple "replace all" in your text editor will do
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,SA=1802,SS=SU+SU,S1=SU+1,SD=SU+SU+1,
     *  ST1=ST+1,ST2=ST+ST,SE=SD+2,SF=SU+SU+SU+SU+2)
C	Sometimes you only need to calculate cross sections, not intensities.
C	In that case, use the line PARAMETER(SAa=1...
C	Advantage?  You will have more memory left, so SU can be higher.
C	Otherwise, for Mueller matrix calculations, do PARAMETER(SAa=SA.. 
C	PARAMETER(SAa=1,SSs=1,SUu=1,SDd=1,S11=1,SEe=1)
	PARAMETER(SAa=SA,SSs=SS,SUu=SU,SDd=SD,S11=S1,SEe=SE)
	REAL*8 EPS,EPS1,EPS2,EPS3,KE,KE1,KE2,KE3,KI,KI1,KI2,KI3,KEX
	REAL*8 CEXT,CSCA,QEXT,QSCA,QABS,ALFA,QBACK,ASIM,CGN1,CGN2
	REAL*8 MR2,DELTA,SEN,COSS,XM,CGN(200),CTE,NH
	REAL*8 AP1(SDd)
	REAL*8 KM,KIM,EM,DPRIMA,PC,PE,PI,H,DN,pc1,pc2,DB(st2+2)
	REAL*8 sum,lam
	real*8 PX1(SSs,SAa),PX2(SSs,SAa),PX3(SSs,SAa),PX4(SSs,SAa)
	real*8 AP2(SDd),AP3(SDd),AP4(SDd),BP1(SDd),BP2(SDd) 
	real*8 AE1(SAa),AE2(SAa),AE3(SAa),AE4(SAa),BE1(SAa),BE2(SAa)
	real*8 A1(SAa),A2(SAa),A3(SAa),A4(SAa),B1(SAa),B2(SAa),ANG(SAa)
	real*8 FALLOS(SAa)
	COMPLEX*16 D3(SUu,SUu,SDd),D4(SUu,SUu,SDd),D5(SUu,SUu,SDd)
	COMPLEX*16 AX1(SUu,SUu,SDd),AX2(SUu,SUu,SDd)
	COMPLEX*16 BX1(SUu,SDd,SDd),BX2(SUu,SDd,SDd)
	COMPLEX*16 G1(SDd),G2(SDd)
	COMPLEX*16 CT,Z1,Z2,Z3,Z4,Z5,TA(2),MR,MR1
	COMPLEX*16 T1(ST2,ST2,ST1),T(ST2,ST2,ST1)
	complex*16 G3(SDd),G4(SDd),G5(SDd)
	INTEGER NOCON1,NOCON2,NOCON,NCERO,STOPE,CFA(5),NC,CF
	INTEGER I,M,N,N1,NPR,NSO,M1,S,NGAUSS,N1MAX,N2MAX,N1MAX1,N1MAX2
	INTEGER IAJ,IBJ,ICJ,IAM,IBM,IAY,CZ,NTOPE,NX,DE,CA,OCP
	EXTERNAL PUJOL,TNATURAL,CGORDANR
	COMMON/CAUNO/N1MAX1
	COMMON/CADOS/N1MAX2
	PI=3.141592653589793238462643D0
	PG=PI/180.0D0
C	mr1= Absolute RI of the core (complex)
C	mr2= Absolute RI of the shell (real, nonabsorbing)
C	Does it mean I cannot use absorbing material for the shell?
C	Not unless the code is modified.  We never used absorbing shell material
C	for our laboratory particles, so I never bothered to do it.
	mr1=(1.2d0,0.01d0)
	mr2=1.05d0
	DO I=1,5
		WRITE (*,*)
	END DO
	WRITE (*,*) '               Welcome to Arturo Quirantes`'
	WRITE (*,*)
	WRITE (*,*) '          ==================================='
	WRITE (*,*) '          ============  L I S A  ============'
	WRITE (*,*) '          ==================================='
	WRITE (*,*)
	WRITE (*,*)
	WRITE (*,*) 'This code calculates light-scattering properties'
	WRITE (*,*) 'for axisymmetric, randomly-oriented coated particles'
	WRITE (*,*)
	WRITE (*,*) 'Please enter Delta accuracy parameter'
	WRITE (*,*) '(see Eq. 21 in Appl.Opt. 32, 4652-4666, 1993)'
	WRITE (*,*) '(1): Delta = 0.01  (higher speed)'
	WRITE (*,*) '(2): Delta = 0.001 (higler accuracy)'
	WRITE (*,*) '(3): User-selected Delta'
	WRITE (*,*) '(4): No delta. NMAX and NGAUSS are entered'
	WRITE (*,*) '     (See reference above)'
10      READ (*,*) DE
C	T is a matrix of dimension (2*n1max)*(2*n1max)
C	N2max is the "dimension" needed to calculate cross sections
C	Since usually N2max<N1MAX, this will speed calculations
C	And NGAUSS is the number of Gauss quadrature points needed
C	for the calculation of the T-matrix elements
C	If DE=4, you can set those elements yourself
C	Otherwise, the PUJOL convergence subroutine will do it for you
C	I personally use option 2, unless I need more accurady (DE=3)
	IF(DE.NE.1.AND.DE.NE.2.AND.DE.NE.3.AND.DE.NE.4) GO TO 10
	IF(DE.EQ.1) DELTA=1.0D-2
	IF(DE.EQ.2) DELTA=1.0D-3
	IF(DE.EQ.3) THEN
20        WRITE (*,*) 'Please enter a value for Delta:'
		READ (*,*) DELTA
		IF(DELTA.LE.0) GO TO 20
	END IF
	IF(DE.EQ.4) THEN
		WRITE (*,*) 'Enter a value for N1MAX'
		READ (*,*) N1MAX
		WRITE (*,*) 'Enter a value for N2MAX'
		READ (*,*) N2MAX
		IF(N2MAX.GT.N1MAX) N2MAX=N1MAX
		WRITE (*,*) 'Enter a value for NGAUSS'
		READ (*,*) NGAUSS
		DELTA=0
	END IF
C	If you can to calculate intensities, polarizations and the like,
C	you have to set the number of angles (NTOPE)
C	in this example, theta=0,10,20...180 degrees	
	NTOPE=7
	NC=1
	DO N=1,NTOPE
		ANG(N)=DFLOAT(N-1)*180.0d0/dfloat(ntope-1)
		ANG(N)=ANG(N)*PG
	END DO
	OPEN(1,FILE='q.dat')
	OPEN(2,FILE='failures.dat')
	OPEN(3,FILE='mm1.dat')
	OPEN(4,FILE='mm2.dat')
	OPEN(11,FILE='mm3.dat')
	WRITE (*,*)
	WRITE (*,*) 'You can input the particle`s dimensions in two ways:'
	WRITE (*,*) '1) Outer and Inner size parameters (K*Rext, K*Rint)'
	WRITE (*,*) '2) Outer size (K*Rext) and core/particle ratio (Rint/Rext)'
	WRITE (*,*) 'What`s your choice? (1/2)'
50	READ (*,*) OCP
	IF(OCP.NE.1.AND.OCP.NE.2) GO TO 50
	WRITE (*,*) 
	WRITE (*,*) 'Now, enter size parameters'
	WRITE (*,*) 'Remember, these are dimensionless size parameters'
	WRITE (*,*) 'Please refer to subroutine TNATURAL for more information'
	WRITE (*,*) 
	WRITE (*,*) 'Enter value for Rext (outer surface)'
	READ (*,*) KEX
	IF (OCP.EQ.1) WRITE (*,*) 'Enter value for K*Rint (inner surface)'
	IF (OCP.EQ.2) WRITE (*,*) 'Enter value for Rint/Rext ratio'
	READ (*,*) KI
	WRITE (*,*) 'Entre the value for the nonsphericity parameter'
	WRITE (*,*) 'Please refer to subroutine TNATURAL for more information'
	READ (*,*) EPS
	NOCON=0
C	Now the T-matrix core (TNATURAL) has to be run twice:
C	once for the inner surface (CA=1), and one for the outer surface (CA=2)
C	If you set NMAX, NGAUSS (DE=4), the code will go right to TNATURAL
C	Otherwise, the convergence subroutine PUJOL will be called
	IF(DE.EQ.4) THEN
		CA=1
		MR=MR1/MR2
		KE=KI*MR2
		IF(OCP.EQ.2) KE=KE*KEX
		NPR=N2MAX
		N2MAX=N1MAX
		CALL TNATURAL(CA,MR,KE,EPS,N1MAX,N2MAX,NGAUSS,CEXT,CSCA,T1)
		N2MAX=NPR
		N1MAX1=N1MAX
		CA=2
		MR=MR2
		KE=KEX
		CALL TNATURAL(CA,MR,KE,EPS,N1MAX,N2MAX,NGAUSS,CEXT,CSCA,T)
	ELSE
C	This is for an automatic convergence (nonzero zeta values)
C	NOCON1,NOCON2 are convergence data (see PUJOL for more details)
		CA=1
		MR=MR1/MR2
		KE=KI*MR2
		IF(OCP.EQ.2) KE=KE*KEX
		CALL PUJOL(CA,MR,KE,EPS,DELTA,CEXT,CSCA,T1,NOCON,N1MAX,N2MAX)
		N1MAX1=N1MAX
		NOCON1=NOCON
		CA=2
		MR=MR2
		KE=KEX
		CALL PUJOL(CA,MR,KE,EPS,DELTA,CEXT,CSCA,T,NOCON,N1MAX,N2MAX)
		NOCON2=NOCON
	END IF
C	The T-matrix calculations are finished
C	Now let's calculate the extinction, scattering and absorption efficiencies
C	Please note that Cext, Csca are NOT cross sections
	QEXT=CEXT*2.0D0/(KEX*KEX)
	QSCA=CSCA*2.0D0/(KEX*KEX)
	QABS=QEXT-QSCA
C	If all you need is Qext, Qsca, Qabs please add the following line:
c	go to 444
C	That way, you'll skip the long and tedious intensity calculations
C	On the other hand, if you want Mueller matrix elements, do nothing more
C	Still here?  Great.
C	Now to the expansion coefficients and orientation averaging
C	I suggest JOSA A 8, 871-882 (from now on, JOSA) for a clearer picture
C	Calculation of Ann'n1 (JOSA, Eqs. 4.18)
C
C	Warning:
	do n=0,n2max+n2max
C	Some people have complained that the last line should read "do n=1..."
C	in order to avoid the db(n=0) case
C	However, note that the integer variable in the array db is n+1, NOT n
C	So even in the n=0 case we have db(0+1)=d(1), and everything's fine
		db(n+1)=dfloat(n+n+1)
	end do
	DO 303 N=1,N2MAX
	DO 302 N1=0,N2MAX+N2MAX
	DO 301 NPR=MAX(1,ABS(N-N1)),MIN(N+N1,N2MAX)
	IAM1=0
	IAM2=MIN(N,NPR)
	IBM=0
	CALL CGORDANR(N,N1,NPR,IAM1,IAM2,IBM,CGN)
		IF(MOD((N+N1+NPR),2).EQ.0) THEN
			AX1(N,NPR,N1+1)=CGN(0-IAM1+1)*(T(N,NPR,1)
     *  			+T(N+N1MAX,NPR+N1MAX,1))*0.5D0
			AX2(N,NPR,N1+1)=CGN(0-IAM1+1)*(T(N,NPR,1)
     *  			-T(N+N1MAX,NPR+N1MAX,1))*0.5D0
		ELSE
			AX1(N,NPR,N1+1)=0.0D0
			AX2(N,NPR,N1+1)=0.0D0
		END IF
	DO 300 M1=1,MIN(N,NPR)
	if(CGN(M1-IAM1+1).ne.0) then
		IF(MOD((N+N1+NPR),2).EQ.0) THEN
			TA(1)=T(N,NPR,M1+1)+T(N+N1MAX,NPR+N1MAX,M1+1)
			TA(2)=T(N,NPR,M1+1)-T(N+N1MAX,NPR+N1MAX,M1+1)
		ELSE
			TA(1)=T(N,NPR+N1MAX,M1+1)+T(N+N1MAX,NPR,M1+1)
			TA(2)=T(N,NPR+N1MAX,M1+1)-T(N+N1MAX,NPR,M1+1)
		END IF
		AX1(N,NPR,N1+1)=AX1(N,NPR,N1+1)+TA(1)*CGN(M1-IAM1+1)
		AX2(N,NPR,N1+1)=AX2(N,NPR,N1+1)+TA(2)*CGN(M1-IAM1+1)
	end if
300	CONTINUE
		CT=2.0D0*(0.0D0,1.0D0)**(NPR-N)/SQRT(db(npr+1))
		AX1(N,NPR,N1+1)=AX1(N,NPR,N1+1)*CT
		AX2(N,NPR,N1+1)=AX2(N,NPR,N1+1)*CT
301	CONTINUE
302	CONTINUE
303	CONTINUE
C	Calculation of Bmnn1 (JOSA, Eqs. 4.17)
	DO 330 N=1,N2MAX
	DO 320 N1=0,N2MAX+N2MAX
	DO M=-N2MAX,N2MAX
		if(BX1(N,N1+1,M+N2MAX+1).ne.0) BX1(N,N1+1,M+N2MAX+1)=0.0D0
		if(BX2(N,N1+1,M+N2MAX+1).ne.0) BX2(N,N1+1,M+N2MAX+1)=0.0D0
		END DO
	DO 310 NPR=MAX(1,ABS(N-N1)),MIN(N+N1,N2MAX)
		IAM1=-N2MAX
		IAM2=N2MAX
		IBM=-1
		CALL CGORDANR(N,NPR,N1,IAM1,IAM2,IBM,CGN)
		H=SQRT(DB(NPR+1)/DB(N1+1))
		DO IAM=IAM1,IAM2
			CGN(IAM-IAM1+1)=CGN(IAM-IAM1+1)*H
			IF(MOD((N-IAM),2).NE.0) CGN(IAM-IAM1+1)=-CGN(IAM-IAM1+1)
		END DO
	DO 309 M=-N2MAX,N2MAX
		if(CGN(M-IAM1+1).ne.0) then
		BX1(N,N1+1,M+N2MAX+1)=BX1(N,N1+1,M+N2MAX+1)+CGN(M-IAM1+1)*AX1(N,NPR,N1+1)
		BX2(N,N1+1,M+N2MAX+1)=BX2(N,N1+1,M+N2MAX+1)+CGN(M-IAM1+1)*AX2(N,NPR,N1+1)
		end if
309	CONTINUE
310	CONTINUE
320	CONTINUE
330	CONTINUE
C	Calculation of Dmnnso (JOSA, Eqs. 4.29-4.33)
	DO 349 N=1,N2MAX
	DO 348 NSO=1,N2MAX
		DO 342 M=-MIN(N,NSO),MIN(N,NSO)
			ax1(N,NSO,M+N2MAX+1)=0.0D0
			ax2(N,NSO,M+N2MAX+1)=0.0D0
			DO 341 N1=ABS(M-1),N2MAX+N2MAX
			if(BX1(N,N1+1,M+N2MAX+1).ne.0.and.BX1(NSO,N1+1,M+N2MAX+1).ne.0) then
				ax1(N,NSO,M+N2MAX+1)=ax1(N,NSO,M+N2MAX+1)+DB(N1+1)*
     *				BX1(N,N1+1,M+N2MAX+1)*CONJG(BX1(NSO,N1+1,M+N2MAX+1))
			end if
			if(BX2(N,N1+1,M+N2MAX+1).ne.0.and.BX2(NSO,N1+1,M+N2MAX+1).ne.0) then
				ax2(N,NSO,M+N2MAX+1)=ax2(N,NSO,M+N2MAX+1)+DB(N1+1)*
     *			BX2(N,N1+1,M+N2MAX+1)*CONJG(BX2(NSO,N1+1,M+N2MAX+1))
			end if
341			CONTINUE
342		CONTINUE
		DO 344 M=MAX(-N,-NSO+2),MIN(N,NSO+2)
			D3(N,NSO,M+N2MAX+1)=0.0D0
			D4(N,NSO,M+N2MAX+1)=0.0D0
			D5(N,NSO,M+N2MAX+1)=0.0D0
		DO 343 N1=ABS(M-1),N2MAX+N2MAX
		if(BX1(N,N1+1,M+N2MAX+1).ne.0.and.BX1(NSO,N1+1,2-M+N2MAX+1).ne.0) then
			D3(N,NSO,M+N2MAX+1)=D3(N,NSO,M+N2MAX+1)+DB(N1+1)*
     *	   	BX1(N,N1+1,M+N2MAX+1)*CONJG(BX1(NSO,N1+1,2-M+N2MAX+1))
		end if
		if(BX2(N,N1+1,M+N2MAX+1).ne.0) then
		if(BX2(NSO,N1+1,2-M+N2MAX+1).ne.0) then
			D4(N,NSO,M+N2MAX+1)=D4(N,NSO,M+N2MAX+1)+DB(N1+1)*
     *	   	BX2(N,N1+1,M+N2MAX+1)*CONJG(BX2(NSO,N1+1,2-M+N2MAX+1))
		if(BX1(NSO,N1+1,2-M+N2MAX+1).ne.0) then
			D5(N,NSO,M+N2MAX+1)=D5(N,NSO,M+N2MAX+1)+DB(N1+1)*
     *	   	BX2(N,N1+1,M+N2MAX+1)*CONJG(BX1(NSO,N1+1,2-M+N2MAX+1))
		end if
		end if
		end if
343		CONTINUE
344		CONTINUE
348	CONTINUE
349	CONTINUE
C	Calculation of the G1-G5 (JOSA Eqs. 4.23 - 4.27)
C	and As,Bs (Id. Eqs. 2.31-2.36)
	DO S=0,N2MAX+N2MAX
		G1(S+1)=0.0D0
		G2(S+1)=0.0D0
		G3(S+1)=0.0D0
		G4(S+1)=0.0D0
		G5(S+1)=0.0D0
		AP1(S+1)=0.0D0
		AP2(S+1)=0.0D0
		AP3(S+1)=0.0D0
		AP4(S+1)=0.0D0
		BP1(S+1)=0.0D0
		BP2(S+1)=0.0D0
	END DO
	DO 410 S=0,N2MAX+N2MAX
	DO 400 N=1,N2MAX
	IF(ABS(N-S).GT.N2MAX) GO TO 400
	DO 390 NSO=MAX(1,ABS(N-S)),MIN(N+S,N2MAX)
	Z1=0.0D0
	Z2=0.0D0
	Z3=0.0D0
	Z4=0.0D0
	Z5=0.0D0
c	To save memory, D1 and D2 are stored at AX1,AX2
	IAM1=0
	IAM2=MIN(N,NSO)
	IBM=0
	CALL CGORDANR(N,S,NSO,IAM1,IAM2,IBM,CGN)
	DO 360 M=-MIN(N,NSO),MIN(N,NSO)
	IF(M.LT.0) THEN
		IF(CGN(-M-IAM1+1).NE.0) THEN
		IF(MOD((N+S-NSO),2).EQ.0) THEN
			Z1=Z1+CGN(-M-IAM1+1)*ax1(N,NSO,M+N2MAX+1)
			Z2=Z2+CGN(-M-IAM1+1)*ax2(N,NSO,M+N2MAX+1)
		ELSE
			Z1=Z1-CGN(-M-IAM1+1)*ax1(N,NSO,M+N2MAX+1)
			Z2=Z2-CGN(-M-IAM1+1)*ax2(N,NSO,M+N2MAX+1)
		END IF
		END IF
	ELSE
		IF(CGN(M-IAM1+1).NE.0) THEN
		IF(CGN1.EQ.0.AND.M.EQ.1) CGN1=CGN(M-IAM1+1)
		Z1=Z1+CGN(M-IAM1+1)*ax1(N,NSO,M+N2MAX+1)
		Z2=Z2+CGN(M-IAM1+1)*ax2(N,NSO,M+N2MAX+1)
	END IF
	END IF
360	CONTINUE
C	D4 and D5 are stored as D1 and D2, to save memory
	IF(S.GE.2) THEN
	IAM1=-MIN(N,NSO+2)
	IAM2=MIN(N,NSO-2)
	IBM=2
	CALL CGORDANR(N,S,NSO,IAM1,IAM2,IBM,CGN)
	DO 380 M=MAX(-N,-NSO+2),MIN(N,NSO+2)
		IF(CGN(-M-IAM1+1).NE.0) THEN
			IF(CGN2.EQ.0.AND.M.EQ.1) CGN2=CGN(-M-IAM1+1)
			Z3=Z3+CGN(-M-IAM1+1)*D3(N,NSO,M+N2MAX+1)
			Z4=Z4+CGN(-M-IAM1+1)*D4(N,NSO,M+N2MAX+1)
			Z5=Z5+CGN(-M-IAM1+1)*D5(N,NSO,M+N2MAX+1)
		END IF
380	CONTINUE
	END IF
	H=SQRT(DB(N+1)/DB(NSO+1))*DB(S+1)/(CSCA+CSCA)
	IF(MOD((N+NSO-S),2).EQ.0) THEN
		G2(S+1)=G2(S+1)+H*Z2*CGN1
		IF(S.GE.2) G4(S+1)=G4(S+1)+H*Z4*CGN2
	ELSE
		G2(S+1)=G2(S+1)-H*Z2*CGN1
		IF(S.GE.2) G4(S+1)=G4(S+1)-H*Z4*CGN2
	END IF
	G1(S+1)=G1(S+1)+H*Z1*CGN1
	IF(S.GE.2) G3(S+1)=G3(S+1)+H*Z3*CGN2
	IF(S.GE.2) G5(S+1)=G5(S+1)-H*Z5*CGN1
	CGN1=0.0D0
	CGN2=0.0D0
390	CONTINUE
400	CONTINUE
	AP1(S+1)=DREAL(G1(S+1)+G2(S+1))
	AP2(S+1)=DREAL(G3(S+1)+G4(S+1))
	AP3(S+1)=DREAL(G3(S+1)-G4(S+1))
	AP4(S+1)=DREAL(G1(S+1)-G2(S+1))
	BP1(S+1)=2.0D0*DREAL(G5(S+1))
	BP2(S+1)=2.0D0*IMAG(G5(S+1))
	XM=ABS(AP1(S+1))
	IF(ABS(AP2(S+1)).GT.XM) XM=ABS(AP2(S+1))
	IF(ABS(AP3(S+1)).GT.XM) XM=ABS(AP3(S+1))
	IF(ABS(AP4(S+1)).GT.XM) XM=ABS(AP4(S+1))
	IF(ABS(BP1(S+1)).GT.XM) XM=ABS(BP1(S+1))
	IF(ABS(BP2(S+1)).GT.XM) XM=ABS(BP2(S+1))
	XM=XM*1000.0D0								 
	IF((DELTA.EQ.0).AND.(XM.LT.1.0D-3)) THEN
		STOPE=S
		GO TO 420
	ELSE IF((DELTA.NE.0).AND.(XM.LT.DELTA)) THEN
		STOPE=S
          GO TO 420
	END IF
410	CONTINUE
		STOPE=N2MAX+N2MAX
420	CONTINUE
C	Calculation of angle-dependent coefficients
	DO N=1,NTOPE
	SEN=SIN(ANG(N))
	COSS=COS(ANG(N))
	PX1(1,N)=COSS
	PX1(2,N)=(3.0D0*COSS*COSS-1.0D0)/2.0D0
	PX2(1,N)=0.0D0
	PX2(2,N)=-SEN*SEN*DSQRT(0.375D0)
	PX3(1,N)=0.0D0
	PX3(2,N)=(1.0D0+COSS)*(1.0D0+COSS)/4.0D0
	PX4(1,N)=0.0D0
	PX4(2,N)=(1.0D0-COSS)*(1.0D0-COSS)/4.0D0
	A1(N)=AP1(1)+AP1(2)*PX1(1,N)+AP1(3)*PX1(2,N)
	A4(N)=AP4(1)+AP4(2)*PX1(1,N)+AP4(3)*PX1(2,N)
	A2(N)=(AP2(3)*(PX3(2,N)+PX4(2,N))+AP3(3)*(PX3(2,N)-PX4(2,N)))/2.0D0
	A3(N)=(AP2(3)*(PX3(2,N)-PX4(2,N))+AP3(3)*(PX3(2,N)+PX4(2,N)))/2.0D0
	B1(N)=BP1(3)*PX2(2,N)
	B2(N)=BP2(3)*PX2(2,N)
	DO S=3,STOPE
		PX1(S,N)=(2.0D0*S-1.0D0)*COSS*PX1(S-1,N)-(S-1.0D0)*PX1(S-2,N)
		PX1(S,N)=PX1(S,N)/S
		PX2(S,N)=(2.0D0*S-1.0D0)*COSS*PX2(S-1,N)
		PX2(S,N)=PX2(S,N)-DSQRT((S-1.0D0)*(S-1.0D0)-4.0D0)*PX2(S-2,N)
		PX2(S,N)=PX2(S,N)/DSQRT(S*S-4.0D0)
		PX3(S,N)=(2.0D0*S-1.0D0)*(S*(S-1.0D0)*COSS-4.0D0)*PX3(S-1,N)
		PX3(S,N)=PX3(S,N)-S*((S-1.0D0)*(S-1.0D0)-4.0D0)*PX3(S-2,N)
		PX3(S,N)=PX3(S,N)/((S-1.0D0)*(S*S-4.0D0))
		PX4(S,N)=(2.0D0*S-1.0D0)*(S*(S-1.0D0)*COSS+4.0D0)*PX4(S-1,N)
		PX4(S,N)=PX4(S,N)-S*((S-1.0D0)*(S-1.0D0)-4.0D0)*PX4(S-2,N)
		PX4(S,N)=PX4(S,N)/((S-1.0D0)*(S*S-4.0D0))
		A1(N)=A1(N)+AP1(S+1)*PX1(S,N)
		A2(N)=A2(N)+(AP2(S+1)*(PX3(S,N)+PX4(S,N))+AP3(S+1)*(PX3(S,N)-
     *   		PX4(S,N)))/2.0D0
		A3(N)=A3(N)+(AP2(S+1)*(PX3(S,N)-PX4(S,N))+AP3(S+1)*(PX3(S,N)+
     *   		PX4(S,N)))/2.0D0
		A4(N)=A4(N)+AP4(S+1)*PX1(S,N)
		B1(N)=B1(N)+BP1(S+1)*PX2(S,N)
		B2(N)=B2(N)+BP2(S+1)*PX2(S,N)
	END DO
C	Now for a few inequalities that have to be fulfilled
C	If any of these fail (to withing an accuracy of DELTA/100),
C	Then there's something wrong the calculation process
C	(Most likely, convergence failures)
C	CFA(1)..CFA(5) will tell you which inequality fails, and for how many angles
	FALLOS(1)=(A1(N)+A2(N))*(A1(N)+A2(N))
	FALLOS(1)=FALLOS(1)-(A3(N)+A4(N))*(A3(N)+A4(N))
	FALLOS(1)=FALLOS(1)-4.0D0*(B1(N)*B1(N)+B2(N)*B2(N))
	FALLOS(1)=FALLOS(1)/(A1(N)*A1(N))
	FALLOS(2)=1.0D0-(A2(N)/A1(N))
	FALLOS(2)=FALLOS(2)-ABS((A3(N)/A1(N))-(A4(N)/A1(N)))
	FALLOS(3)=1.0D0-(B1(N)/A1(N))
	FALLOS(3)=FALLOS(3)-ABS((A2(N)/A1(N))-(B1(N)/A1(N)))
	FALLOS(4)=1.0D0+(B1(N)/A1(N))
	FALLOS(4)=FALLOS(4)-ABS((A2(N)/A1(N))+(B1(N)/A1(N)))
C	DPRIMA is the accuracy parameter allowed to the inequalities
C	If you want to change it, here's your chance
	DPRIMA=DELTA/100
	IF(DELTA.EQ.0) DPRIMA=1.0D-6
	IF(FALLOS(1).LT.(-DPRIMA)) CFA(1)=CFA(1)+1
	IF(FALLOS(2).LT.(-DPRIMA)) CFA(2)=CFA(2)+1
	IF(FALLOS(3).LT.(-DPRIMA)) CFA(3)=CFA(3)+1
	IF(FALLOS(4).LT.(-DPRIMA)) CFA(4)=CFA(4)+1
	END DO
	CF=CFA(1)+CFA(2)+CFA(3)+CFA(4)+CFA(5)
C	
C	If the code has come this far, congratulations, you've made it
C	The output parameters are A1(n)..A4(n),B1(n),B2(n)
C	That means that, for an angle ANG(n), the Mueller matrix is:
C
C		(A1 B1  0   0)				(F11 F12  0   0 )
C		(B1 A2  0   0)	or, in different	(F12 F22  0   0 )
C		(0  0   A3 B2)	notation ....		( 0   0  F33 F34)
C		(0  0  -B2 A4)				( 0   0 -F34 F44)
C
C	This is the case for axisymmetric, randomly-oriented particles
C	See e.g. Bohren and Huffman's book
C	Of course, if you chose NOT to calculate any Mueller element...
444	continue
C
C	And here are the final results
C	First, the efficiencies
	WRITE (*,*) 'The results are:'
	WRITE (*,*) 'Qext=',Qext
	WRITE (*,*) 'Qsca=',Qsca
	WRITE (*,*) 'Qabs=',qabs
	WRITE (*,*) 'NOCON convergence data:',nocon1,nocon2
	WRITE (1,*) qext,qsca,qabs
C	Then, the warning in case the inequalities have failed.
	IF(CF.NE.0) then
		WRITE (*,*)
		WRITE (*,*) 'W A R N I N G'
		WRITE (*,*) 'Some checking inequalities have failed'
		WRITE (*,*) 'Please see the efficiencies.dat file for details'
		write (*,*)
		write (2,*) cfa(1),fallos(1)
		write (2,*) cfa(2),fallos(2)
		write (2,*) cfa(3),fallos(3)
		write (2,*) cfa(4),fallos(4)
	END IF
1000	CONTINUE
	WRITE (*,*) 'Want Mueller matrix data stored in file? (1=yes, 0=no)'
1010	READ (*,*) NCERO
	IF(NCERO.NE.1.AND.NCERO.NE.0) GO TO 1010
	IF (NCERO.EQ.1) THEN
	DO S=1,NTOPE
		ANG(S)=ANG(S)/PG			 
		WRITE (3,*) S,REAL(ANG(S)),A1(S),A2(S)
		WRITE (4,*) S,REAL(ANG(S)),A3(S),A4(S)
		WRITE (11,*) S,REAL(ANG(S)),B1(S),B2(S)
	END DO
	WRITE (*,*) 'Mueller matrix elements stored on files mm1.dat, mm2.dat and and mm3.dat'
	END IF
1500	WRITE (*,*)
	WRITE (*,*) 'End of the program'
	END
C
C
C
C
C	SUBROUTINE FOR CONVERGENCE IN N1MAX,N2MAX,NGAUSS
	SUBROUTINE PUJOL(CA,MR,KE,EPS,DELTA,CEXT,CSCA,T,NOCON,
     *  N1MAX,N2MAX)      
	INTEGER SU,SUP,ST,ST1,ST2
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,ST1=ST+1,ST2=ST+ST)
	REAL*8 CEXT,CSCA,K1,K2,K3,K4,K5,K6,DELTA,C1(2),C2(2),KE,EPS
	COMPLEX*16 MR,CCX,CBX,T(ST2,ST2,ST1)
	INTEGER*4 NOCON
	INTEGER NMAX,N1MAX,N2MAX,NGAUSS,ML,CA,N1MAX2,GP
	EXTERNAL TNATURAL
	COMMON/CADOS/N1MAX2
	NOCON=0
C       ========== Convergence for N1max ==========
C	We need a starting value for NMAX
C	The first logical step is Wiscombe's criterion x+4x^0.3333+2
C	However, this usually means a great degree of accuracy
C	So, it's better to lower a bit, by using ADHOC
C	Try your own values.  You can try to raise it, if needed.
	ADHOC=-10
	if(ca.eq.2) ADHOC=-15
10	NMAX=INT(KE+4.0D0*((KE)**(1.0D0/3.0D0))+2)
	NMAX=NMAX+ADHOC
C	Of course, we do want NMAX to be positive, hmm?
	IF(NMAX.LT.1) NMAX=1
C	And we also want to have NMAX(core)<NMAX(shell)
	IF(CA.EQ.2.AND.NMAX.LT.N1MAX2) NMAX=N1MAX2
C	GP is the step size for Gauss integration
C	The number of integration points is set as NMAX*GP
C	GP=4 is a good value, but you might want to raise it for
C	highly nonspherical particles.
	GP=4
	NGAUSS=NMAX*GP
C	For convergence checking, we only need m=0
	ML=0
	CALL TNATURAL(CA,MR,KE,EPS,NMAX,ML,NGAUSS,CEXT,CSCA,T)
	C1(2)=CEXT
	C2(2)=CSCA
	K3=1000.0D0
	K4=1000.0D0
	K1=100.0D0
	K2=100.0D0
30	NMAX=NMAX+1
C	If NMAX is larger than we can afford, the we must stop
	IF(NMAX.GT.ST) THEN
		WRITE (*,*) 'N1MAX IS HIGHER THAT THE MAXIMUM POSSIBLE VALUE OF',SU
		STOP
	END IF
	NGAUSS=NMAX*GP
	CALL TNATURAL(CA,MR,KE,EPS,NMAX,ML,NGAUSS,CEXT,CSCA,T)
	C1(1)=C1(2)
	C2(1)=C2(2)
	C1(2)=CEXT
	C2(2)=CSCA
	K1=ABS(1-C1(1)/C1(2))
	K2=ABS(1-C2(1)/C2(2))
C	K1, K2 are the differences between Cext, Csca for Nmax and for Nmax-1
C	(To be precise, Cext and Csca are only calculated using m=0)
C	If K1,K2<0.1*Delta, then our convergence goal has been achieved
C	Otherwise, we make NMAX+1, and continue
	IF(MAX(K1,K2).GT.DELTA/10) GO TO 30
40	CONTINUE
	N1MAX=NMAX
C       ========== Convergence for N2max ==========
C	We have calculated the T matrix for m=0 to within the desired accuracy
C	However, many times we do not need too use all their elements
C	So we will get N2MAX, and will not use T-matrix elements with a larger index
C	If CA=1 (inner surface), it's better to use the complete matrix: N2max=N1max
	IF(CA.EQ.1) THEN
		N2MAX=N1MAX
		GO TO 60
	END IF
C	For the outer surface, let's see how many T-matrix elements are needed
C	We will calculate Cext, Csca for N2max and for N1max
C	When those sections converge to within 0.1*Delta, we will stop
C	A good starting point will be half N1max
	N2MAX=MAX(INT(N1MAX/2),1)
C	Let's now calculate Cext,Csca
C	We will again suppose m=0
	CBX=0.0D0
	CCX=(0.0D0,0.0D0)
	DO I=1,N2MAX
		DO J=1,N2MAX
			CBX=CBX+ABS(T(I,J,1))*ABS(T(I,J,1))
			CBX=CBX+ABS(T(I+NMAX,J+NMAX,1))*ABS(T(I+NMAX,J+NMAX,1))
		END DO
		CCX=CCX+T(I,I,1)+T(I+NMAX,I+NMAX,1)
	END DO
	C1(2)=-DREAL(CCX)
	C2(2)=CBX
50	N2MAX=N2MAX+1
C	If N2MAX=N1MAX, it means we need all T-matrix elements; well...
	IF(N2MAX.EQ.N1MAX) GO TO 60
	CBX=0.0D0
	CCX=(0.0D0,0.0D0)
	DO I=1,N2MAX
		DO J=1,N2MAX
			CBX=CBX+ABS(T(I,J,1))*ABS(T(I,J,1))
			CBX=CBX+ABS(T(I+NMAX,J+NMAX,1))*ABS(T(I+NMAX,J+NMAX,1))
		END DO
		CCX=CCX+T(I,I,1)+T(I+NMAX,I+NMAX,1)
	END DO
	C1(2)=-DREAL(CCX)
	C2(2)=CBX
	K1=ABS(1-C1(2)/CEXT)
	K2=ABS(1-C2(2)/CSCA)
C	Just like before, if we get convergente better than 0.1*Delta, success!
	IF(MAX(K1,K2).GT.DELTA/10) GO TO 50
60	CONTINUE
C       ========== Convergence in NGAUSS ==========
C	And now, let's see how many Gauss quadrature points we need
C	We will start with GP*NMAX
	NMAX=N1MAX
	ML=0
	NGAUSS=NMAX*GP
	CALL TNATURAL(CA,MR,KE,EPS,NMAX,ML,NGAUSS,CEXT,CSCA,T)
	K3=1000.0D0
	K4=1000.0D0
	K1=100.0D0
	K2=100.0D0
70	C1(1)=CEXT
	C2(1)=CSCA
C	And now, let's increase NGAUSS by GP
	NGAUSS=NGAUSS+GP
	CALL TNATURAL(CA,MR,KE,EPS,NMAX,ML,NGAUSS,CEXT,CSCA,T)
C	And now for a little difference
C	Here I demand that the differences in Cext,Csca decreases
C	for two consecutive values of Ngauss:
C	That is, Cext(ngauss)<Cext(ngauss-1)<Cext(ngauss-2)
C	The reason is, K1 and K2 can increase before going down again
C	So I let them go up on the condition that they go down again
C	If they don't, it might mean the convergence cannot be achieved
C	And further increases in Ngauss might only worsen results
C	In the past, I did the same check in N1MAX, but since it cannot
C	go higher than SU, it would have to stop anyway.
	K5=K3
	K6=K4
	K3=K1
	K4=K2
	K1=ABS(1-C1(1)/CEXT)
	K2=ABS(1-C2(1)/CSCA)
	IF((MAX(K1,K2).GT.MAX(K3,K4)).AND.(MAX(K3,K4).GT.MAX(K5,K6))) THEN
		NOCON=NOCON+10000000
		NGAUSS=NGAUSS-GP-GP
		GO TO 80
	END IF
	IF(MAX(K1,K2).GT.DELTA/10) GO TO 70
80	CONTINUE
	CALL TNATURAL(CA,MR,KE,EPS,N1MAX,N2MAX,NGAUSS,CEXT,CSCA,T)
	NOCON=INT(NGAUSS+1000*N2MAX+100000*N1MAX)+NOCON
C	Please take a look at NOCON.  It's a reminder of the T-matrix size
C	It has the form abbccddd, where bb=N1max, cc=N2max, dd=Ngauss
C	a=1 if there is no convergence in Ngauss, 0 otherwise
	IF(CA.EQ.1) N1MAX2=N1MAX
	RETURN
	END
C
C
C
C
C	SUBROUTINE FOR CALCULATION OF THE NATURAL T-MATRIX
C	(Natural=Z axis along the symmetry axis)
C	This is the core of the apple, so to speak
C	I've made it along the lines of Barber and Hill: "Light Scattering
C	by Particles: Computational Methods" (Barber-Hill from now on)
C	With some changes e.g. the use of Wigner instead of Legendre functions
C	So don't be surprised if the equations slightly differ from that book's
C	In case of doubt, read my lips: it works!
	SUBROUTINE TNATURAL(CA,MR,KE,EPS,NMAX,ML,NGAUSS,CEXT,CSCA,T)
	INTEGER SU,SUP,ST,ST1,ST2,ST3
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,ST1=ST+1,ST2=ST+ST,ST3=ST2+1)
	DOUBLE PRECISION SENTH,COSTH,SENT2,I1,DERI,THETA,WTSEN
	DOUBLE PRECISION DESEN,KE,EPS,RX(200),RW(200),WIG(ST1),WIG2(ST1)
	DOUBLE PRECISION PI,KR,IPJ,JA(ST1),JZ(ST1),GM(ST1)
	DOUBLE PRECISION DEL(ST1),D(ST1),DP(ST1),ZB,ZE,ZFX,ZFY,GM1
	DOUBLE PRECISION KAYUDA,CEXT,CSCA,DELT
	DOUBLE COMPLEX X1,ZAX,ZDX,IA(ST1),IZ(ST1),IB(ST1),KA(ST1),CB(ST1)
	DOUBLE COMPLEX A(ST2,ST2),B(ST2,ST2),AA(ST2,ST2),BB(ST2,ST2)
	DOUBLE COMPLEX HA(ST1),BC(ST1),BC2(ST1)
	DOUBLE COMPLEX ZA,ZC,ZD,ZF,ZR
	DOUBLE COMPLEX GM2,MR,MD,MKR
	DOUBLE COMPLEX T(ST2,ST2,ST1),T1(ST2,ST2,ST1),CCX
	INTEGER REL(ST1),REL2(ST1)
	INTEGER M,I,J,IREL,JREL,TBUCLE,N1MAX1,k,nq
	INTEGER N,MTOPE,R,NMAX,ML,NGAUSS,SIM,CA,NN,NK
	EXTERNAL BESSEL,BESSEL2,HANKEL,WIGNER,GAUSS
	EXTERNAL LULABY
	COMMON/CAUNO/N1MAX1
	PI=3.141592653589793238462643D0
C       CA=1: core.     CA=2: coating.
C	if m1=m2, we have a homgeneous scatterer, then T1=0 and over for CA=1
C	Same if the core has zero size
	if(ca.eq.1.and.(mr.eq.1.or.ke.eq.0)) then
		do m=0,ml
			MTOPE=MAX(1,M)
			R=M+1
			N=NMAX-MTOPE+1
			DO I=1,N+N
			DO J=1,N+N
				IF(T1(I,J,R).NE.0) T1(I,J,R)=(0.0D0,0.0D0)
			END DO
			END DO
		end do
		cext=1.0d0
		csca=1.0d0
		go to 120
	end if
C       Gauss quadrature subroutine.  Rx,Rw: abscissas and weights
	CALL GAUSS(NGAUSS,RX,RW)
	MD=MR*MR
C	If the particle has a plane of symmetry perpendicular to the symmetry
C	axis, calculations get simpler.  Half of the T-matrix elements need not be
C	computed (they equal zero), and the other half need surface integration
C	only from 0 to pi/2
C	For the cases of plane-symmetric particles, SIM=1, otherwise SIM=0
	SIM=1
C	Now, let's calculate n and n*(n+1) in hight-precision
	DO I=1,NMAX+1
		D(I)=DFLOAT(I)
		DP(I)=DFLOAT(I)*DFLOAT(I+1)
		if(I.NE.0) GM(I)=SQRT((2.0D0*D(I)+1.0D0)/DP(I))
	END DO
C	Loop for M (azimuthal parameter)
C	That means, the T-matrix decomposes into m T-submatrices
C	That happens for axisymmetric particles.
	DO 40 M=0,ML
	MTOPE=MAX(1,M)
	R=M+1
	N=NMAX-MTOPE+1
	DO I=MTOPE,NMAX
		REL(I)=I-MTOPE+1
	END DO
C       Loop for TBUCLE, (angle)
	DO 30 TBUCLE=1,NGAUSS
C	As said before, if SIM=1 we need not integrate from pi/2 to pi
C	In that case, the value of the A,B matrices is halved
C	But we will calculate B*A^-1, so there's really no need to
C	multiply by 2 for each A,B, matrix element
C	Maybe there's a more elegant means to halve Ngauss prior to the
C	GAUSS subroutine, but this will also do
		IF(SIM.EQ.1.AND.TBUCLE.GT.(NGAUSS/2)) GO TO 30
		THETA=RX(TBUCLE)
		SENTH=SIN(THETA)
		SENT2=SENTH*SENTH
		COSTH=COS(THETA)
		WTSEN=RW(TBUCLE)*SENTH
C	Subroutine to calculate Wigner functions
C	They don't explode for hight n,m values like Legendre functions
C	WIG(N+1)=dm0_n
		CALL WIGNER (SENTH,COSTH,M,NMAX,WIG)
C	Now we have to specify the particle surface: r(theta) and derivative
C	This is easy, as long as the particle is axisymmetric
C	I've used a spheroid, and of course you can put your own shape
C	You can even specify different surfaces for core and coating
C	e.g. an inner spheroid inside a cylinder.
C	All you have to add is an "IF CA=1 ... else ..."
C	The limit is your imagination.
C	Of course, I suggest you make sure that the core fits INSIDE the particle!
C
C	Spheroidal particle with axes a,b (a=axis of revolution)
C	KE=Equivalent-volume size parameter = (a*b*b)^0.33333, eps=b/a
		SIM=1
		KR=KE*EPS**(1.0D0/3.0D0)
		KR=KR/SQRT(SENTH*SENTH+(COSTH*EPS)**2)
		DERI=SENTH*COSTH*(EPS*EPS-1.0D0)*KR*KR*KR
		DERI=DERI/(KE*KE*EPS**(2.0D0/3.0D0))
C	Calculation of Hankel and Bessel function
C	HA(n+1)=hn(kr), BC(n+1)=jn(mkr)
		MKR=MR*KR
		CALL HANKEL(KR,NMAX,HA)
		CALL BESSEL(MKR,NMAX,BC)
C	And in the outer surface, we also need Hankel functions for mkr
		IF(CA.EQ.2) THEN
			CALL BESSEL2(MKR,NMAX,BC2)
			DO I=1,NMAX+1
				BC2(I)=(0.0D0,1.0D0)*BC2(I)+BC(I)
			END DO
		END IF
		DESEN=DERI*SENTH
C	Now for another variation
C	In Barber-Hill (pp. 90-91), Legendre polynomials are used in pairs:
C	P(i)*P(j), P(i-1)*P(j), P(i)*P(j-1), P(i-1)*P(j-1)
C	First, they were substituted for Wigner functions
C	But there's a problem: when RX nears unity, P's are also close to one
C	So, ZE (see below) was calculated as a difference of nearly-unity terms
C	Therefore increasing the chance of roundoff errors
C	I've since introduced a Delta function as:
C	Del(i)=Cos(theta)*Del(i-1)-i*Sen(theta)*Sen(theta)*Wigner(i)
C	(that's for m=0; when m<>0, its slightly different
C	That way, calculations are simplified and can be done with a bit more
C	accuracy
	IF(M.EQ.0) THEN
		DEL(1)=-SENT2
		DEL(2)=-(SENT2+SENT2+SENT2)*COSTH
		DO I=3,NMAX
			DEL(I)=DEL(I-1)*COSTH-D(I)*SENT2*WIG(I)
		END DO
	ELSE
		DO I=MTOPE,NMAX
			DEL(I)=D(I)*COSTH*WIG(I+1)-WIG(I)*SQRT((D(I)-D(M))*(D(I)+D(M)))
		END DO
	END IF
C	And now, let's calculate A, B matrix elements
C	I=row, J=column
C	See that A,B =0 when (i or j are less than m) and m>1
C	So the loop goes from MTOPE to NMAX, MTOPE being the minimum value
C	of i,j so that we only calculate nonzero matrix elements
		DO I=MTOPE,NMAX
			IZ(I)=HA(I+1)*WTSEN
			IA(I)=HA(I)/HA(I+1)
			JZ(I)=DREAL(IZ(I))
			JA(I)=DREAL(HA(I))/DREAL(HA(I+1))
			IB(I)=MR*BC(I)/BC(I+1)
			KA(I)=JA(I)-IA(I)
			WIG2(I)=WIG(I+1)*DP(I)*DESEN	
			IF(CA.EQ.2) CB(I)=MR*BC2(I)/BC2(I+1)
		END DO
		DO 20 I=MTOPE,NMAX
 		DO 10 J=MTOPE,NMAX 
			IPJ=D(I)*D(J)/KR
			I1=WIG(I+1)*WIG(J+1)*DESEN
			GM1=0.5D0*GM(I)*GM(J)/SENT2
			IF(MOD(M,2).NE.0) GM1=-GM1
			GM2=(0.0D0,1.0D0)*GM1
C	Let's now use the following sub-matrix description for A,B
C	First quadrant (I) is the (nmax*nmax) left upper side of the matrix
C	Second quadrant (II) is the (nmax*nmax) right upper side of the matrix
C	Third (III) and fourth (IV) are the same for the left and right lower side
C
C	II and III calculations.
C	For plane-symmetric particles and i+j=even, matrix elements are zero
		IF((SIM.EQ.1.AND.MOD((I+J),2).NE.0).OR.SIM.EQ.0) THEN
C     And, if M=0, II and III quadrant elements equal zero
		IF(M.NE.0) THEN
			ZA=KR+IA(I)*(KR*IB(J)-D(J))-D(I)*IB(J)+IPJ
			ZB=DEL(I)*WIG(J+1)+DEL(J)*WIG(I+1)
			ZB=ZB*KR
			ZC=I1*(DP(I)*IB(J)+DP(J)*IA(I)-IPJ*(D(I)+D(J)+2.0D0))
			ZR=ZC+KA(I)*DP(J)*I1
			X1=D(M)*GM1*BC(J+1)
C	======== III quadrant calculations ============
	A(N+REL(I),REL(J))=A(N+REL(I),REL(J))-IZ(I)*X1*(ZA*ZB+ZC)
 			ZAX=ZA+KA(I)*(KR*IB(J)-D(J))
	B(N+REL(I),REL(J))=B(N+REL(I),REL(J))-JZ(I)*X1*(ZAX*ZB+ZR)
C	======== II quadrant calculations  ============
			ZA=ZA+KR*(MD-1.0D0)
	A(REL(I),N+REL(J))=A(REL(I),N+REL(J))-IZ(I)*X1*(ZA*ZB+ZC)/MR
			ZAX=ZA+KA(I)*(KR*IB(J)-D(J))
	B(REL(I),N+REL(J))=B(REL(I),N+REL(J))-JZ(I)*X1*(ZAX*ZB+ZR)/MR
C       =============================================================
C       =============================================================
C       Additional matrices for outer surface (CA=2); first part
		IF(CA.EQ.2) THEN
			ZA=KR*(1.0D0+IA(I)*CB(J))-D(J)*IA(I)-D(I)*CB(J)+IPJ
			ZC=ZC+(CB(J)-IB(J))*DP(I)*I1
			ZR=ZC+KA(I)*DP(J)*I1		
			X1=D(M)*GM1*BC2(J+1)
	AA(N+REL(I),REL(J))=AA(N+REL(I),REL(J))-IZ(I)*X1*(ZA*ZB+ZC)
			ZAX=ZA+KA(I)*(KR*CB(J)-D(J))		
	BB(N+REL(I),REL(J))=BB(N+REL(I),REL(J))-JZ(I)*X1*(ZAX*ZB+ZR)
			ZA=KR*(MD+IA(I)*CB(J))-D(J)*IA(I)-D(I)*CB(J)+IPJ
	AA(REL(I),N+REL(J))=AA(REL(I),N+REL(J))-IZ(I)*X1*(ZA*ZB+ZC)/MR
			ZAX=ZA+KA(I)*(KR*CB(J)-D(J))		
	BB(REL(I),N+REL(J))=BB(REL(I),N+REL(J))-JZ(I)*X1*(ZAX*ZB+ZR)/MR
		END IF
C       =============================================================
C       =============================================================
		END IF
		END IF                  
C	I and IV quadrants calculations
C	For plane-symmetric particles and i+j=odd, matrix elements are zero
		IF((SIM.EQ.1.AND.MOD((I+J),2).EQ.0).OR.SIM.NE.1) THEN
			ZD=KR*(IB(J)-MD*IA(I))+MD*D(I)-D(J)
			ZE=DEL(I)*DEL(J)
			IF(M.NE.0) ZE=ZE+(DP(M)-D(M))*WIG(I+1)*WIG(J+1)
			ZE=ZE*KR				
			ZFX=WIG2(J)*DEL(I)
			ZFY=WIG2(I)*DEL(J)
			X1=GM2*BC(J+1)
C       ======== IV quadrant calculations ==============
			ZF=ZFX-MD*ZFY
	A(N+REL(I),N+REL(J))=A(N+REL(I),N+REL(J))+IZ(I)*X1*(ZD*ZE+ZF)/MR
			ZDX=ZD-(JA(I)-IA(I))*KR*MD
	B(N+REL(I),N+REL(J))=B(N+REL(I),N+REL(J))+JZ(I)*X1*(ZDX*ZE+ZF)/MR
C       ======== I quadrant calculations  ==============
			ZD=KR*(IB(J)-IA(I))+D(I)-D(J)
			ZF=ZFX-ZFY
	A(REL(I),REL(J))=A(REL(I),REL(J))+IZ(I)*X1*(ZD*ZE+ZF)
			ZDX=ZD-(JA(I)-IA(I))*KR
	B(REL(I),REL(J))=B(REL(I),REL(J))+JZ(I)*X1*(ZDX*ZE+ZF)
C       =============================================================
C       =============================================================
C       Additional matrices for outer surface (CA=2); second part
	IF(CA.EQ.2) THEN
			ZD=KR*(CB(J)-MD*IA(I))+MD*D(I)-D(J)
			ZF=ZFX-MD*ZFY
			X1=GM2*BC2(J+1)
	AA(N+REL(I),N+REL(J))=AA(N+REL(I),N+REL(J))+IZ(I)*X1*(ZD*ZE+ZF)/MR
			ZDX=ZD-(JA(I)-IA(I))*KR*MD
	BB(N+REL(I),N+REL(J))=BB(N+REL(I),N+REL(J))+JZ(I)*X1*(ZDX*ZE+ZF)/MR
			ZD=KR*(CB(J)-IA(I))+D(I)-D(J)
			ZF=ZFX-ZFY
	AA(REL(I),REL(J))=AA(REL(I),REL(J))+IZ(I)*X1*(ZD*ZE+ZF)
			ZDX=ZD-(JA(I)-IA(I))*KR
	BB(REL(I),REL(J))=BB(REL(I),REL(J))+JZ(I)*X1*(ZDX*ZE+ZF)
		END IF
C       =============================================================
C       =============================================================
		END IF
C	End of the J (column) loop
10	CONTINUE         
C     	End of the I (row) loop
20	CONTINUE
C	End of the TBUCLE (surface angle) loop
30	CONTINUE
C	Now some matrix multiplication for CA=2
	IF(CA.EQ.2) THEN
	NN=N1MAX1-MTOPE+1
	NK=MIN(N,NN)
	IF(NK.EQ.0) GO TO 45
	DO 43 I=1,NN
	DO 42 J=1,NN
	DO 41 K=1,NK
	IF((SIM.EQ.1.AND.MOD((I+J),2).EQ.0).OR.SIM.NE.1) THEN
		A(I,J)=A(I,J)+AA(I,K)*T1(K,J,R)+AA(I,K+N)*T1(K+NN,J,R)
		B(I,J)=B(I,J)+BB(I,K)*T1(K,J,R)+BB(I,K+N)*T1(K+NN,J,R)
		A(I+N,J+N)=A(I+N,J+N)+AA(I+N,K)*T1(K,J+NN,R)+AA(I+N,K+N)*T1(K+NN,J+NN,R)
		B(I+N,J+N)=B(I+N,J+N)+BB(I+N,K)*T1(K,J+NN,R)+BB(I+N,K+N)*T1(K+NN,J+NN,R)
	END IF
	IF(M.NE.0) THEN
	IF((SIM.EQ.1.AND.MOD((I+J),2).NE.0).OR.SIM.NE.1) THEN
		A(I,J+N)=A(I,J+N)+AA(I,K)*T1(K,J+NN,R)+AA(I,K+N)*T1(K+NN,J+NN,R)
		B(I,J+N)=B(I,J+N)+BB(I,K)*T1(K,J+NN,R)+BB(I,K+N)*T1(K+NN,J+NN,R)        
		A(I+N,J)=A(I+N,J)+AA(I+N,K)*T1(K,J,R)+AA(I+N,K+N)*T1(K+NN,J,R)
		B(I+N,J)=B(I+N,J)+BB(I+N,K)*T1(K,J,R)+BB(I+N,K+N)*T1(K+NN,J,R)
	END IF
	END IF
41	CONTINUE
42	CONTINUE
43	CONTINUE
	END IF
45	CONTINUE
C	Enough with multiplications
C	Now we have to invert+multiply:  T = -B*A^(-1)
C	If CA=1, T is the T1 matrix; otherwise, if CA=2, T is the final T-matrix
C	The LULABY matrix does the product B*A^(-1), storing it in B
C	nq is the size of the A, B, matrix to invert
	nq=nmax-max(1,m)+1
	nq=nq+nq
	CALL LULABY(A,B,NQ,ca,ml,m)
	DO J=1,NMAX+NMAX
	DO I=1,NMAX+NMAX
		IF(CA.EQ.1.AND.T1(I,J,R).NE.0) T1(I,J,R)=(0.0D0,0.0D0)
		IF(T(I,J,R).NE.0) T(I,J,R)=(0.0D0,0.0D0)
	END DO
	END DO
C	And now we "decompress" the T (or T1) matrix
	DO I=1,N
	DO J=1,N
		IREL=I+MTOPE-1
		JREL=J+MTOPE-1
		T(IREL,JREL,R)=-B(I,J)
		T(IREL+NMAX,JREL,R)=-B(I+N,J)
		T(IREL,JREL+NMAX,R)=-B(I,J+N)
		T(IREL+NMAX,JREL+NMAX,R)=-B(I+N,J+N)
		IF(CA.EQ.1) THEN
			T1(I,J,R)=-B(I,J)
			T1(I+N,J,R)=-B(I+N,J)
			T1(I,J+N,R)=-B(I,J+N)
			T1(I+N,J+N,R)=-B(I+N,J+N)
		END IF
		IF (A(I,J).NE.0) A(I,J)=(0.0D0,0.0D0)
		IF (B(I,J).NE.0) B(I,J)=(0.0D0,0.0D0)
		IF (A(I+N,J).NE.0) A(I+N,J)=(0.0D0,0.0D0)
		IF (B(I+N,J).NE.0) B(I+N,J)=(0.0D0,0.0D0)
		IF (A(I,J+N).NE.0) A(I,J+N)=(0.0D0,0.0D0)
		IF (B(I,J+N).NE.0) B(I,J+N)=(0.0D0,0.0D0)
		IF (A(I+N,J+N).NE.0) A(I+N,J+N)=(0.0D0,0.0D0)
		IF (B(I+N,J+N).NE.0) B(I+N,J+N)=(0.0D0,0.0D0)
		IF(CA.EQ.2) THEN
			IF (AA(I,J).NE.0) AA(I,J)=(0.0D0,0.0D0)
			IF (BB(I,J).NE.0) BB(I,J)=(0.0D0,0.0D0)
			IF (AA(I+N,J).NE.0) AA(I+N,J)=(0.0D0,0.0D0)
			IF (BB(I+N,J).NE.0) BB(I+N,J)=(0.0D0,0.0D0)
			IF (AA(I,J+N).NE.0) AA(I,J+N)=(0.0D0,0.0D0)
			IF (BB(I,J+N).NE.0) BB(I,J+N)=(0.0D0,0.0D0)
			IF (AA(I+N,J+N).NE.0) AA(I+N,J+N)=(0.0D0,0.0D0)
			IF (BB(I+N,J+N).NE.0) BB(I+N,J+N)=(0.0D0,0.0D0)
		END IF
	END DO
	END DO
C	Enf of the azimuth (m) loop
40	CONTINUE
C	Calculation of extinction and scattering "sections"
	CSCA=0.0D0
	CEXT=0.0D0
	CCX=(0.0D0,0.0D0)
100	CONTINUE
	DO M=0,ML
		MTOPE=MAX(M,1)-1
		do i=1,nmax
		do j=1,nmax
			IF (M.EQ.0) THEN
				DELT=1.0D0
			ELSE
				DELT=2.0D0
			END IF
		KAYUDA=(ABS(T(I,J,M+1)))*(ABS(T(I,J,M+1)))
		KAYUDA=KAYUDA+(ABS(T(I+NMAX,J,M+1)))*(ABS(T(I+NMAX,J,M+1)))
		KAYUDA=KAYUDA+(ABS(T(I,J+NMAX,M+1)))*(ABS(T(I,J+NMAX,M+1)))
		KAYUDA=KAYUDA+(ABS(T(I+NMAX,J+NMAX,M+1)))*(ABS(T(I+NMAX,
     *		J+NMAX,M+1)))
		CSCA=CSCA+DELT*KAYUDA
		KAYUDA=0.0D0    
		END DO
		END DO
	END DO
	ccx=(0.0d0,0.0d0)
	DO I=1,NMAX
		CCX=CCX+(T(I,I,1)+T(I+NMAX,I+NMAX,1))
	END DO
	DO M=1,ML
		DO I=M,NMAX
			CCX=CCX+2.0D0*(T(I,I,M+1)+T(I+NMAX,I+NMAX,M+1))
		END DO
	END DO
110	CONTINUE
	CEXT=-DREAL(CCX)
C	If you want Cross Sections, multiply CEXT, CSCA by 2*pi/(k*k)
C	where k = 2*pi/wavelength
C	Or, to get Qext, Qsca, multiply Cext, Csca by 2/(X*X)
C	where X is the dimensionless size parameter
C	You can do it in the main routine.  Here we're through
120	RETURN
	END
C
C
C
C
C
C	Subroutine to calculate the complex-argument Bessel j(n) functions
C	by downward recurrence.	BESSEL(n)=BC(n+1)
	SUBROUTINE BESSEL(MKR,NMAX,BC)
	INTEGER SU,SUP,ST,ST1,SX
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,ST1=ST+1,SX=ST1+100)
	DOUBLE COMPLEX QB,BC(ST1),MKR,XC(SX),XZ,BZ
	DOUBLE PRECISION N2
	INTEGER T,N1
	N2=MAX(NMAX,INT(ABS(MKR)))
	N1=INT(N2+4.0D0*N2**0.333333+8.0d0*SQRT(N2)+5)
	XC(N1)=(0.0D0,0.0D0)
	XC(N1-1)=(1.0D-35,0.0D0)
	DO T=N1-2,1,-1
		XC(T)=DFLOAT(T+T+3)*XC(T+1)/MKR-XC(T+2)
	END DO
	XZ=3.0D0*XC(1)/MKR-XC(2)
	BZ=SIN(MKR)/MKR
	QB=BZ/XZ
	BC(1)=BZ 
	DO T=2,NMAX+1
		BC(T)=XC(T-1)*QB
	END DO
	CONTINUE
	RETURN
	END
C
C
C
C
C
C	Subroutine to calculate the complex-argument Bessel y(n) functions
C	by upward recurrence.  BESSEL2(n)=BC2(n+1)
	SUBROUTINE BESSEL2(MKR,NMAX,BC2)
	INTEGER SU,SUP,ST,ST1,SX
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,ST1=ST+1)
	DOUBLE COMPLEX BC2(ST1),MKR
	INTEGER T
	BC2(1)=-COS(MKR)/MKR
	BC2(2)=(BC2(1)-SIN(MKR))/MKR
	DO T=2,NMAX
		BC2(T+1)=DFLOAT(T+T-1)*BC2(T)/MKR-BC2(T-1)
	END DO
	CONTINUE
	RETURN
	END     
C	
C
C
C
C
C	Subroutine to calculate the real-argument Hankel function
C	by upward and downward recurrence.  HANKEL(n)=HA(n+1)
	SUBROUTINE HANKEL(KR,NMAX,HA)
	INTEGER SU,SUP,ST,ST1,SX
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,ST1=ST+1,SX=ST1+100)
	DOUBLE PRECISION QA,JJ(SX),YR(ST1),KR,JZ,N2
	DOUBLE COMPLEX HA(ST1)
	INTEGER T,N1
	N2=MAX(NMAX,INT(ABS(KR)))
	N1=INT(N2+4.0D0*N2**0.333333+8.0d0*SQRT(N2)+5)
	JJ(N1)=0.0D0
	JJ(N1-1)=1.0D-35
	DO T=N1-2,1,-1
		JJ(T)=DFLOAT(T+T+3)*JJ(T+1)/KR-JJ(T+2)
	END DO
	JZ=3.0D0*JJ(1)/KR-JJ(2)
	QA=(SIN(KR)/KR)/JZ
	YR(1)=-COS(KR)/KR
	YR(2)=(YR(1)-SIN(KR))/KR
	HA(1)=DCMPLX((SIN(KR)/KR),YR(1))
	HA(2)=DCMPLX(JJ(1)*QA,YR(2))
	DO T=2,NMAX
		YR(T+1)=DFLOAT(T+T-1)*YR(T)/KR-YR(T-1)
		HA(T+1)=DCMPLX(JJ(T)*QA,YR(T+1))
	END DO
	CONTINUE
	RETURN
	END
C
C
C
C
C
C	Subroutine to calculate the Wigner functions
C	WIG(n+1)=dm0,n
	SUBROUTINE WIGNER (SENTH,COSTH,M,NMAX,WIG)
	INTEGER SU,SUP,ST,ST1
	PARAMETER (SU=40,SUP=0,ST=SU+SUP,ST1=ST+1)
	INTEGER TI
	DOUBLE PRECISION SENTH,COSTH,WIG(ST1)
C	M=0
	IF(M.EQ.0) THEN
		WIG(1)=1.0D0
		WIG(2)=COSTH
		DO TI=2,NMAX
			WIG(TI+1)=(DFLOAT(TI+TI)-1.0D0)*COSTH*WIG(TI)
			WIG(TI+1)=WIG(TI+1)-(DFLOAT(TI)-1.0D0)*WIG(TI-1)
			WIG(TI+1)=WIG(TI+1)/(DFLOAT(TI))
		END DO
		RETURN
C	M<>0
	ELSE
		WIG(M+1)=1.0D0
		IF(MOD(M,2).NE.0) WIG(M+1)=-WIG(M+1)
		DO TI=0,M-1
			WIG(TI+1)=0.0D0
			WIG(M+1)=WIG(M+1)*SENTH*DSQRT(DFLOAT(TI+1)*DFLOAT(M+TI+1))
			WIG(M+1)=WIG(M+1)/(2.0D0*DFLOAT(TI+1))
		END DO
		DO TI=M+1,NMAX
			WIG(TI+1)=DFLOAT(TI+TI-1)*COSTH*WIG(TI)
			WIG(TI+1)=WIG(TI+1)-DSQRT(DFLOAT((TI-1)*(TI-1)-M*M))*WIG(TI-1)
			WIG(TI+1)=WIG(TI+1)/DSQRT(DFLOAT(TI*TI-M*M))
		END DO
	RETURN
	END IF
	CONTINUE
	END
C
C
C
C
C
C	=========================================================
C	Subroutine for LU-based inversion and multiplication 
C	of a square nq*nq matrix
C	The final result (B) is equal to the product (A^-1)*B
C	Since we need to do B*(A^-1) we transpose the A, B matrices:
C	And then we will transpose the final matrix:
C	(At^-1)*Bt = (B*(A^-1))t = -Tt
C	This is a taylored adaptation of the ZGESV (=ZGETRF+ZGETRS)
C	subroutine from the LAPACK package
	SUBROUTINE LULABY(A,B,NQ,ca,ml,m)
	INTEGER SU,SUP,ST2,ST3
	PARAMETER (SU=40,SUP=0,ST2=SU+SU+SUP+SUP)
	DOUBLE COMPLEX a(ST2,ST2),b(ST2,ST2),uno,cero,ztemp
	DOUBLE PRECISION SMAX		  
	INTEGER IPIV(ST2),NQ,J,K,LL,IMAX,JP,INFO,IP,ca,ml,m
	uno=(1.0d+0,0.0d+0)
	cero=(0.0d+0,0.0d+0)
	info=0
C	First of all, transpose the A, B, matrices
	do j=1,nq
	do i=1,j-1
		ztemp=a(i,j)
		a(i,j)=a(j,i)
		a(j,i)=ztemp
		ztemp=b(i,j)
		b(i,j)=b(j,i)
		b(j,i)=ztemp
	end do
	end do
C	ZGETRF routine to make an LU decomposition
	DO 100 J=1,NQ
C	Search for pivoting and checking for singularity
C	Sub-subroutine IZAMAX to get the value of i such that a(i,j)=max for a given j
		nt=nq-j+1
		imax=1
		if(nt.eq.1) go to 20
                smax=cdabs(a(j,j))
		do 10 i=2,nt
			if(abs(a(j+i-1,j)).le.smax) go to 10
			imax=i
			smax=cdabs(a(j+i-1,j))
10		continue
20		continue
		nt=nq-j
C	End of sub-subroutine IZAMAX
		jp=j-1+imax
		ipiv(j)=jp
		if(a(jp,j).ne.cero) then
C	Sub-subroutine ZSWAP for row interchange
			if(jp.ne.j) then
				do i=1,nq
					ztemp=a(j,i)
					a(j,i)=a(jp,i)
					a(jp,i)=ztemp
				end do
C	End of sub-subroutine ZSWAP
			end if
C	Sub-subroutine ZSCAL to rescale: divide by por a(j,j)
			if(j.lt.nq) then
				ztemp=uno/a(j,j)
				do i=1,nt
					a(j+i,j)=ztemp*a(j+i,j)
				end do
C	End of sub-subroutine ZSCAL
			end if
C	Now if info=1 that means that U(l,l)=0 and invertion of A is not feasible
		else if(info.eq.0) then
			info=j
			write (*,*) 'MATRIZ U SINGULAR EN info=',J
			stop
		end if			 
		if(j.lt.nq) then
C	Sub-subroutine ZGERU para substitute A by A*x*y' (x,y vectors)
C	And therefore update the 'trailing submatrix'
		do k=1,nt
			if(a(j,j+k).ne.cero) then
				ztemp=-uno*a(j,j+k)
				do ll=1,nt
					a(j+ll,j+k)=a(j+ll,j+k)+a(j+ll,j)*ztemp
				end do
			end if
		end do 
C	End of sub-subroutine ZGERU
		end if
C	End of the J loop J
100	continue
C	End of subroutine ZGETRF
C
C	A is now an LU decomposition of a rowwise permutation of A
C	Diagonal unity elements belong tl T.  That is, Lii=1
C
C	Resolving to obtain X = A^-1 * B
c	The procedure is to solve L*X=B, and then U*X=B
C	Sub-subroutine ZLASWP for row interchange in B	
	do i=1,nq
		ip=ipiv(i)
		if(ip.ne.i) then
			do k=1,nq
				ztemp=b(i,k)
				b(i,k)=b(ip,k)
				b(ip,k)=ztemp
			end do
		end if
c	End of subroutine ZLASWP
	end do	
c	Subroutine ZTRSM to solve L*X=B.  The result is overwritten in B
	do j=1,nq
		do k=1,nq
			if(b(k,j).ne.cero) then
				do i=k+1,nq
					b(i,j)=b(i,j)-b(k,j)*a(i,k)
				end do
			end if
		end do
	end do
c	Subroutine ZTRSM to solve U*X=B.  The final result is overwritten in B
	do j=1,nq
		do k=nq,1,-1
			if(b(k,j).ne.cero) then
				b(k,j)=b(k,j)/a(k,k)
				do i=1,k-1
					b(i,j)=b(i,j)-b(k,j)*a(i,k)
				end do
			end if
		end do
	end do
C	End of subroutine ZTRSM
C	And now, let's transpose again
C	The final result, B, equals the product B*A^(-1)
	do j=1,nq
	do i=1,j-1
		ztemp=b(i,j)
		b(i,j)=b(j,i)
		b(j,i)=ztemp
	end do
	end do	 
C	End of subroutine LULABY
	return
	end
C
C
C
C
C	Subroutine to calculate weights and abscissas in Gauss integration.
C	Reference: "Light Scattering by Particles: Computational Methods",
C	P.W. Barber y S.C. Hill, pp. 172-173
	SUBROUTINE GAUSS(NGAUSS,RX,RW)
	DOUBLE PRECISION RX(200),RW(200),RN1(200),RN2(200)
	DOUBLE PRECISION CN,CNN1,CON1,CON2,XI,X,PM1,PM2,RN,AUX,DER1P
	DOUBLE PRECISION DER2P,APPFCT,B,BISQ,BFROOT,RATIO,PROD
	DOUBLE PRECISION CONST,TOL,C1,C2,C3,C4,PI,P,PMP1,PP
	INTEGER NDIV2,NP1,NM1,NM2,IN,NGAUSS,K
	DATA PI/3.141592653589793238462643D0/
	DATA CONST/0.148678816357D0/
	DATA TOL/1.0D-15/
	DATA C1,C2/0.125D0,-0.0807291666D0/
	DATA C3,C4/0.2460286458D0,-1.824438767D0/
C	Warning: if TOL is lower than your computer's roundoff error,
C	your code might fall into an endless loop!
	IF(NGAUSS.EQ.1) THEN
		RX(1)=0.577350269189626D0
		RW(1)=1.0D0
		RETURN
	END IF
	CN=DFLOAT(NGAUSS)
	NDIV2=NGAUSS/2
	NP1=NGAUSS+1
	CNN1=CN*(CN+1.0D0)
	APPFCT=1.0D0/SQRT((CN+0.5D0)**2+CONST)
	CON1=0.5D0*PI
	CON2=0.5D0*PI
	do in=2,ngauss
		rn=dfloat(in)
		rn2(in)=(rn-1.0d0)/rn
		rn1(in)=rn2(in)+1.0d0
	end do
	DO 1030 K=1,NDIV2
		B=(DFLOAT(K)-0.25D0)*PI
		BISQ=1.0D0/(B*B)
C	Bisq is a first approximation to Rx
		BFROOT=B*(1.0D0+BISQ*(C1+BISQ*(C2+BISQ*(C3+C4*BISQ))))
		XI=COS(APPFCT*BFROOT)
1010		X=XI
		PM2=1.0D0
		PM1=X
		DO 1020 IN=2,NGAUSS
			RN=DFLOAT(IN)
			P=((2.0D0*RN-1.0D0)*X*PM1-(RN-1.0D0)*PM2)/RN
			PM2=PM1
			PM1=P
1020		CONTINUE
		PM1=PM2
		AUX=1.0D0/(1.0D0-X*X)
		DER1P=CN*(PM1-X*P)*AUX
		DER2P=((X+X)*DER1P-CNN1*P)*AUX
		RATIO=P/DER1P
		XI=X-RATIO*(1.0D0+RATIO*DER2P/(2.0D0*DER1P))
		IF(ABS(XI-X).GT.TOL) GO TO 1010
C	Now a bit of modification here to calculate Pn with higher accuracy
C	The above calculates PM1=Pn-1(X), not Pn-1(XI)
C	If X and XI differ, it might enlarge roundoff errors
C	So just to play in the safe side...
		X=XI
		PM2=1.0D0
		PM1=X
		PMP1=1.0D0
		DO IN=2,NGAUSS
			RN=DFLOAT(IN)
			P=RN1(IN)*X*PM1-RN2(IN)*PM2
			PM2=PM1
			PM1=P
			PP=X*PMP1+RN*PM2
			PMP1=PP
		END DO
C	End of the fine-tuning modification
C	The weights are usually calculated as wi=2*(1-xi^2)/(n*Pn-1)^2
C	However, for small xi values, the relative error of the weight is higher
C	Another way to calculate it is wi=2/[(1-xi^2)(Pn')^2]
C	This works better at the endpoints
		RX(K)=-XI
		RW(K)=2.0d0/((1.0D0-XI*XI)*PMP1*PMP1)
		RX(NP1-K)=-RX(K)
		RW(NP1-K)=RW(K)
1030		CONTINUE
		IF(MOD(NGAUSS,2).NE.0) THEN
		RX(NDIV2+1)=0.0D0
		NM1=NGAUSS-1
		NM2=NGAUSS-2
		PROD=CN
		DO 1040 K=1,NM2,2
			PROD=PROD*DFLOAT(NM1-K)/DFLOAT(NGAUSS-K)
1040		CONTINUE
		RW(NDIV2+1)=2.0D0/(PROD*PROD)
		END IF
	DO 1050 K=1,NGAUSS
		RX(K)=CON1*RX(K)+CON2
		RW(K)=CON1*RW(K)
1050	CONTINUE
	RETURN
	END
C
C
C
C
C
C	Subroutine for Clebsch-Gordan calculations
C	Since we need a huge lot of calculations, this sobroutine
C	does not calculate just one.  Instead, a set of CG are
C	calculated by using recurrende relations.
C	Specifically, we will calculate C-G(aj,bj,cj;am,bm,cm)
C	For given values of aj,bj,cj,bm and all possible values of am
C	(possible means all am values so that C-G is nonzero, of course)
C	It has been checked to Messiah "Quantum mechanics" Appendix C, Eq. 14a
C	The sum is thumbs-up to within a 10^-12 factor when N2MAX=30
C	Reference: Schulten y Gordon, J. Math. Phys. 16, 1961-1970 (1975)
	SUBROUTINE CGORDANR(IAJ,IBJ,ICJ,IAM1,IAM2,IBM,CGN)
	real*8 cc1,cc2,dd1,dd2,cgn(200),CX(200),C1,SC,CMAX,CMAX2
	INTEGER IAJ,IBJ,ICJ,IBM,IAM1,IAM2,IAM,IAMEDIO
	INTEGER DN1,KP,KP2,IAMA,IAMB,DNTOPE,IAMIN,IAMAX,TIC
C	First, let's check that icj and ibm falls within limits
	if(icj.lt.abs(iaj-ibj).or.icj.gt.(iaj+ibj)) go to 100
	if(abs(ibm).gt.ibj) go to 200
C	IAMIN, IAMAX are the lowest and highest possible values for am
	IAMIN=-MIN(IAJ,ICJ+IBM)
	IAMAX=MIN(IAJ,ICJ-IBM)
C	IAMIN, IAMAX might be such that CG=0
C	If so happens, let's jump up till CG<>0
	KP=MAX(0,IAMIN-IAM1)
	KP2=MAX(0,IAM2-IAMAX)
	DNTOPE=10
	TIC=0
C	Calculation of every possible CG value
	IAMA=IAMIN
	IAMB=IAMAX
	DN1=IAMAX-IAMIN
	CX(1)=1.0D0
	CMAX=CX(1)
	SC=CX(1)*CX(1)
	if(dn1.eq.0) go	to 100
	dd2=dfloat(iaj*(iaj+1)-ibj*(ibj+1)+icj*(icj+1))
	cc2=dfloat((iaj-iama+1)*(iaj+iama))
	cc2=cc2*dfloat((icj-iama-ibm+1)*(icj+iama+ibm))
	cc2=-sqrt(cc2)
C	If there are not many CG to calculate, upward recurrence alone will do
	IF(DN1.GE.DNTOPE) IAMB=IAMIN+INT(DN1/2)-1
C	Upward recurrence
	iam=iama
5	iam=iam+1
		icm=iam+ibm
		if(abs(icm).gt.icj) go to 8
		cc1=cc2
		cc2=dfloat((iaj-iam+1)*(iaj+iam))
		cc2=cc2*dfloat((icj-icm+1)*(icj+icm))
		cc2=-sqrt(cc2)
c	C-G(iaj,iaj,icj/iam,iam,icm)=0 if icj=odd, so...
		if(iaj.eq.ibj.and.iam.eq.ibm.and.mod(icj,2).ne.0) then
			CX(IAM-IAMIN+1)=0.0d0
			go to 5
		end if
		dd1=dd2-2.0d0*dfloat((iam-1)*(icm-1))
		CX(IAM-IAMIN+1)=-CX(IAM-IAMIN)*DD1/CC2
		if(dn1.eq.1) then
			SC=SC+CX(IAM-IAMIN+1)*CX(IAM-IAMIN+1)
			go to 100
		end if
		if(iam.eq.(iama+1)) go to 7
		CX(IAM-IAMIN+1)=CX(IAM-IAMIN+1)-CX(IAM-IAMIN-1)*CC1/CC2
7		SC=SC+CX(IAM-IAMIN+1)*CX(IAM-IAMIN+1)
		CMAX2=ABS(CX(IAM-IAMIN+1))
		IF(CMAX2.GT.CMAX) CMAX=CMAX2
C	Now let's see if we need go on with upward recurrence
		if(iam.lt.iamb) go to 5
		if(tic.eq.1) go to 8
c	CX(middle) is used for re-scaling.  If equals zero, then disaster!
c	So let us just check and, if zero, go for the next one.
c	C-G(iaj,iaj,icj/iam,iam,icm)=0 if icj=odd, so...
	if(iaj.eq.ibj.and.iam.eq.ibm.and.mod(icj,2).ne.0) go to 5
c	Due to round-off errors, we have to be a bit more tolerant: X<10^-14
c	(X being the smallest/largest C-G ratio in our series)
c	However, for high values of iaj,ibj... many C-G are very small, and it
c	becomes difficult to know whether C-G is zero because of roundoff errors
c	or whether that is its real value.
c	My solution: use the next C-G in the series by "jumping" just once
	if(dn1.ge.dntope.and.iam.ge.iamb.and.(cmax2/cmax).le.(1.0D-14)) then
		tic=1	
		go to 5
	end if
8	continue
	if(dn1.lt.dntope) go to 100
	iamedio=iam
	C1=CX(IAMEDIO-IAMIN+1)
	SC=SC-C1*C1
C	Downward recurrence
	IAMA=IAMEDIO
	iamb=iamax
	CX(IAMB-IAMIN+1)=1.0D0
	cc1=dfloat((iaj-iamb)*(iaj+iamb+1))
	cc1=cc1*dfloat((icj-iamb-ibm)*(icj+iamb+ibm+1))
	cc1=-sqrt(cc1)
	do iam=iamb-1,iama,-1
		icm=iam+ibm
		if(icm.gt.icj) go to 100
		cc2=cc1
		cc1=dfloat((iaj-iam)*(iaj+iam+1))
		cc1=cc1*dfloat((icj-icm)*(icj+icm+1))
		cc1=-sqrt(cc1)		
c	C-G(iaj,iaj,icj/iam,iam,icm)=0 if icj=odd, so...
	if(iaj.eq.ibj.and.iam.eq.ibm.and.mod(icj,2).ne.0) then
		CX(IAM-IAMIN+1)=0.0d0
		go to 18
	end if
		if(-icm.gt.icj) go to 18
		dd1=dd2-2.0d0*dfloat((iam+1)*(icm+1))
		CX(IAM-IAMIN+1)=-CX(IAM-IAMIN+2)*DD1/CC1
		if(iam.eq.(iamb-1)) go to 18
		CX(IAM-IAMIN+1)=CX(IAM-IAMIN+1)-CX(IAM-IAMIN+3)*CC2/CC1
18	end do
	C1=C1/CX(IAMEDIO-IAMIN+1)
C	Now, rescaling
	do iam=IAMEDIO,IAMAX
		cx(iam-iamin+1)=cx(iam-iamin+1)*c1
		sc=sc+cx(iam-iamin+1)*cx(iam-iamin+1)
	end do
C	Here's where the combined (upward+downward) recurrence ends
100	CONTINUE
	SC=(dfloat(2*ICJ+1)/dfloat(2*IBJ+1))/SC
        SC=SQRT(SC)
C	Now, let's set the sign for SC
	if(cx(iamax-iamin+1).ne.abs(cx(iamax-iamin+1))) SC=-SC
	if(mod((iaj+iamax),2).ne.0) SC=-SC
	if(kp.ne.0) then
		do iam=0,kp-1
			cgn(iam+1)=0.0d0				
		end do
	end if
	if(kp2.ne.0) then
		do iam=iam2-kp2+1,iam2
			cgn(iam-iam1+1)=0.0d0
		end do
	end if
	do iam=iam1+kp,iam2-kp2
		cgn(iam-iam1+1)=cx(iam-iamin+1)*SC
	end do
200	CONTINUE
	RETURN
	END