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: 17 december 2020
C	
C	This is BART, a computer code to calculate light-scattering properties
C	for coated spheres, following the Aden-Kerker formulation
C	It is adapted for polydisperse particle systems
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 for profitable purposes
CC
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
	DOUBLE COMPLEX M1,M2,SU1(2000),SU2(2000),zzc
	DOUBLE PRECISION D1BUCLE,D2BUCLE,RGR,TH(2000),ITEO(2000)
	DOUBLE PRECISION IEE,ICE,NR,LINF,LSUP,CSCA,CEXT,PE,PC,DE1,DE2,CBACK
	DOUBLE PRECISION PG,LONDA,X,Y,GU(2000),Z,XOK,ZOK,SOK,SIGMIN,SIGMAX
	DOUBLE PRECISION SIGVAR,SIGMA,C,XY,DELT,LK,QEXT,QSCA,FU(2003),QBACK
	DOUBLE PRECISION D1MIN,D1MAX,D1VAR,D2MIN,D2MAX,D2VAR,PT,CIN,zzz,zz2
	INTEGER NAS,NCERO,P,PP,TX,OPCION,DD,K,T,CAL,NPI,CONV,CHECKINT(2003)
	COMMON/AK/X,Z,XY,D2BUCLE,M1,M2,SIGMA,P,PP,TH,OPCION,CONV    
	EXTERNAL IRENE

	PG=3.141592653589793238462643D0
	RGR=PG/180.0D0
C	Londa=vacuum wavelength, in nanometers
	LONDA=550.0D0
C	NR=vacuum index of refraction
	NR=1.0D0
C	m1=core index of refraction
	m1=(1.489d0,0.078d0)
	M1=M1/NR
C	m2=shell index of refraction
	m2=(1.02d0,0.0d0)
	M2=M2/NR
C	PT=constant to pass from dimensioness size parameter to radius (nm)
C	If you want to input size parameter, just PT=1
C	Otherwise, PT=2.0D0*PG*NR/LONDA
	PT=1.0D0
	open(1,file='q.dat')
	OPEN(3,FILE='s11.dat')
	open(4,file='s33.dat')
	open(5,file='s12.dat')
	open(6,file='s34.dat')
	WRITE (*,*) '               Welcome to Arturo Quirantes`'
        WRITE (*,*)
	WRITE (*,*) '          ==================================='
	WRITE (*,*) '          ============  B A R T  ============'
	WRITE (*,*) '          ==================================='
	WRITE (*,*)
	WRITE (*,*)
	WRITE (*,*) 'This code calculates light-scattering properties'
	WRITE (*,*) 'for polydisperse, spherical, coated particles'
	WRITE (*,*) 'It can also be used for monodisperse particles'
	WRITE (*,*)
	WRITE (*,*)
	WRITE (*,*) 'Do you want to calculate 1) Cross Sections, 2) Intensities ?'
50      READ (*,*) CAL
	IF (CAL.NE.1.AND.CAL.NE.2) GO TO 50
	WRITE (*,*) 'From now on X=core size, Z=shell size, Y=X+Z=particle size, Q=Z/Y'
	WRITE (*,*) 'One of these parameters will be affected by polydispersity'
	WRITE (*,*)
	WRITE (*,*) 'Please choose what parameter you want to keep MONODISPERSE'
	WRITE (*,*) '1) Monodisperse core X; polydisperse coating Z'
	WRITE (*,*) '2) Monodisperse coating Z; polydisperse core X' 
	WRITE (*,*) '3) Monodisperse relative coating Q; polydisperse particle size Y'
	WRITE (*,*)
100     READ (*,*) OPCION
	IF (OPCION.NE.1.AND.OPCION.NE.2.AND.OPCION.NE.3) GO TO 100
	WRITE (1,*) 'NR,WAVELENGTH:',REAL(NR),REAL(LONDA)
	WRITE (1,*) 'M1:',M1
	WRITE (1,*) 'M2:',M2
	IF (OPCION.EQ.1) WRITE (1,*) 'OPTION 1: CONSTANT CORE'
	IF (OPCION.EQ.2) WRITE (1,*) 'OPTION 2: CONSTANT ABSOLUTE SHELL'
	IF (OPCION.EQ.3) WRITE (1,*) 'OPTION 3: CONSTANT RELATIVE SHELL'
C       Option 1: constant core
C       Option 2: constant shell
C       Option 3: constant core/particle 
	IF (OPCION.EQ.1.OR.OPCION.EQ.2) THEN
		WRITE (*,*) 'Please enter core size parameter value X'
		READ (*,*) D1BUCLE
		WRITE (*,*) 'Please enter shell size parameter value Z'
		READ (*,*) D2BUCLE
	ELSE
		WRITE (*,*) 'Please enter particle size parameter value Y'
		READ (*,*) D1BUCLE
		WRITE (*,*) 'Please enter relative coating Q'
		READ (*,*) D2BUCLE
	END IF
	WRITE (*,*)
	WRITE (*,*) 'Now for the polydisperse part'
	WRITE (*,*) 'We will be using a log-normal size distribution:'
	WRITE (*,*)
	WRITE (*,*) 'p(d)= Constant*(1/d)*exp(-(ln(d)-ln(dm))^2 / (2*sigma)^2 )'
	WRITE (*,*)
	WRITE (*,*) 'Please enter a sigma value (if monodisperse, make sigma=0)'
	WRITE (*,*) 'Warning: if sigma is large, this can take a looooong time'
	READ (*,*) SIGMA
	WRITE (*,*) 'Now choose accuracy parameter'
	WRITE (*,*) 'Case 1: Delta=10-6 for all cases'
	WRITE (*,*) 'Case 2: Delta=10-6 for number of integration points < 1025'
	WRITE (*,*) '        Delta=10-4 for number of integration points < 1025'
160     READ (*,*) CONV
	IF(CONV.NE.1.AND.CONV.NE.2) GO TO 160
170     CONTINUE
	WRITE (*,*)
	WRITE (*,*) 'Data input finished, calculation in process'
	WRITE (*,*)
	PP=3
C       Option cal=1: cross sectionss. PP=2: CEXT, CSCA
	IF(CAL.EQ.1) THEN
		P=0
	ELSE
C       Option cal=2: intensity values
C       P=Number of angles
		P=181	
		DO T=1,P
			TH(T)=DFLOAT(T-1)*180.0d0/dfloat(P-1)
			TH(T)=TH(T)*RGR
		END DO
	END IF
	PC=1.0D6
c       Constant core or shell
	IF (OPCION.EQ.1.OR.OPCION.EQ.2) THEN
		X=D1BUCLE*PT
		Z=D2BUCLE*PT
		Y=X+Z
	ELSE
		Y=PT*D1BUCLE
		X=Y*D2BUCLE
		Z=Y-X
	END IF
C       Parameter affected by polydispersity:
C       Option 1: Z     Option 2: X     Option 3: Y
	IF (OPCION.EQ.1) XY=Z
	IF (OPCION.EQ.2) XY=X
	IF (OPCION.EQ.3) XY=Y
C       If there is no polydispersity.....
	IF(SIGMA.EQ.0) THEN
		LINF=XY
		LSUP=LINF
		CALL IRENE(LINF,LSUP,DE1,DE2,P,PP,NPI,CONV,FU,su1,su2)
                NPI=INT(Y+4.0D0*Y**0.33333+2)
	ELSE
C       Else if there is polydispersity, there's a couple things to do
C	First, we will cut the integration limits from (zero-infinity)
C	To LINF,LSUP
		C=10
		DELT=1.0D-6
		LK=EXP(SQRT(2.0D0*C*LOG(10.0D0))*SIGMA)
		LINF=XY/LK
		LSUP=XY*LK
		DE1=DELT
		IF(CONV.EQ.2) THEN
			DE2=DE1*100.0D0
		ELSE
			DE2=DE1
		END IF
c	Next, let's use the integration subroutine
		CALL IRENE(LINF,LSUP,DE1,DE2,P,PP,NPI,CONV,FU,su1,su2)
	END IF
	IF(CAL.EQ.1) THEN
		CEXT=0.0D0
		CSCA=0.0D0
	ELSE
		DO K=1,P
			ITEO(K)=0.0D0
		END DO
	END IF
	IF(SIGMA.EQ.0) THEN
		CIN=1.0D0
	ELSE
		CIN=1.0D0/(SQRT(2.0D0*PG)*SIGMA)       
	END IF
	IF(P.NE.0) THEN
		DO T=1,P
			ITEO(T)=FU(T)
		END DO
	END IF
C	Here are the cross efficiency values
	QEXT=FU(P+1)*CIN*2.0d0/(y*y)
	QSCA=FU(P+2)*CIN*2.0d0/(y*y)
	WRITE (1,*) 'QEXT',QEXT
	WRITE (1,*) 'QSCA',QSCA
	WRITE (1,*) 'QABS',QEXT-QSCA
800     CONTINUE
C	And here are the values of the nonzero Müller matrix elements
C	For notation, see e.g. Bohren and Huffman's "Absorption and Scattering
C	of Light by Small Particles", page 112
	do t=1,p
C	This is the S11=S22 matrix element
		zz2=(abs(su2(t))*abs(su2(t))+abs(su1(t))*abs(su1(t)))/2.0d0
		write (3,*) REAL(TH(t)/RGR),zz2						 
C	This is the S12=S21 matrix element
		zzz=(abs(su2(t))*abs(su2(t))-abs(su1(t))*abs(su1(t)))/2.0d0
		write (5,*) REAL(TH(t)/RGR),zzz
C	This is the S33=S44 matrix element
		zzc=(conjg(su2(t))*su1(t)+su2(t)*conjg(su1(t)))/2.0d0
		write (4,*) REAL(TH(t)/RGR),real(zzc)
C	This is the  S34=-S43 matrix element
		zzc=(su1(t)*conjg(su2(t))-su2(t)*conjg(su1(t)))
		zzc=zzc*(0.0d0,1.0d0)
		write (6,*) REAL(TH(t)/RGR),real(zzc)
C	The rest of the Müller matrix elements are zero, so we're finished
	end do	
2200    WRITE (*,*)                          
	IF(CAL.EQ.1) GO TO 2500
2500    CONTINUE
	WRITE (*,*) 'Calculation process finished'
	WRITE (*,*)
	WRITE (*,*) 'BART quits now. Go away man!'
	END
C
C
C
C
C       ROMBERG INTEGRATION SUBROUTINE
C       Source: "Numerical analysys.  The mathematics of scientific calculation"
C       (David Kincaid y Ward Cheney).
C	Now I use other subroutines (like Gauss), but this is one I tried
C	In order to improve timing.
C	I still don´t know if it´s faster or not.  Anywa, it works.
	SUBROUTINE IRENE(LINF,LSUP,DE1,DE2,P,PP,NPI,CONV,FU,su1,su2)
	double complex su1(2000),su2(2000)
	DOUBLE PRECISION AY(2000),R1(2000,100),R2(2000,100),POM,PR1,PR2
	DOUBLE PRECISION H,XSUBI,FU(2003),DELTA,LINF,LSUP,DE1,DE2
	INTEGER RN,RM,PO,T,P,PP,NPI,TP,I,CHECKINT(2000),CONV
	EXTERNAL ADENKERKER
	RN=0
	PO=1
	DO T=1,P+PP
		R1(T,1)=0.0D0
	END DO
	H=LSUP-LINF
	DO TP=1,2
		IF(TP.EQ.1) XSUBI=LSUP
		IF(TP.EQ.2) XSUBI=LINF
999             CALL ADENKERKER(XSUBI,FU,SU1,SU2)
		IF(LINF.EQ.LSUP) GO TO 1200
		DO T=1,P+PP
			R1(T,1)=R1(T,1)+FU(T)
		END DO
	END DO
	DO T=1,P+PP
		R1(T,1)=R1(T,1)*H/2.0D0
	END DO
1000    RN=RN+1
	PO=PO*2
	IF(RN.EQ.1) PO=1
	H=H/2.0D0
	IF(CONV.EQ.2.AND.RN.GT.11) THEN 
		DELTA=DE2
	ELSE
		DELTA=DE1
	END IF
	DO T=1,P+PP
		AY(T)=0.0D0
	END DO
	DO I=1,PO
		XSUBI=LINF+(2.0D0*DFLOAT(I)-1.0D0)*H
1001            CALL ADENKERKER(XSUBI,FU,SU1,SU2)
		DO T=1,P+PP
			AY(T)=AY(T)+FU(T)
		END DO
	END DO
	DO T=1,P+PP
		R2(T,1)=AY(T)*H+R1(T,1)/2.0D0
	END DO
	POM=1.0D0
	DO RM=1,RN
		POM=POM*4.0D0
		DO T=1,P+PP
			R2(T,RM+1)=(R2(T,RM)-R1(T,RM))/(POM-1.0D0)
			R2(T,RM+1)=R2(T,RM+1)+R2(T,RM)
		END DO 
	END DO
	DO RM=RN,1,-1
		PR2=0.0D0
		DO T=1,PP
		PR1=ABS(1.0D0-R1(P+T,RM)/R2(P+T,RM))
			IF(PR1.GT.PR2) PR2=PR1
		END DO
		IF(PR2.LT.DELTA) THEN
			IF(CHECKINT(RM).EQ.1) THEN
			NPI=1000*(PO+PO+1)+(RM-1)
			DO T=1,RM
				CHECKINT(RM)=0
				END DO
			GO TO 1100
			ELSE
				CHECKINT(RM)=1
			END IF
		ELSE
			CHECKINT(RM)=0
			DO T=1,P+PP
				R1(T,RM)=R2(T,RM)
			END DO
		END IF
	END DO
	DO T=1,P+PP
		R1(T,RN+1)=R2(T,RN+1)
	END DO
	GO TO 1000
1100	CONTINUE
	DO T=1,P+PP
		FU(T)=R2(T,RM)
	END DO
	RETURN
1200		CONTINUE
	END
C
C
C
C
C
C	This is the core subroutine for Aden-Kerker calculations.
C	I followed the equations as in Bohren and Huffman´s
C	"Absorption and scattering of light by small particles." etc, etc.
	SUBROUTINE ADENKERKER(XSUBI,FU,SU1,SU2)
	DOUBLE COMPLEX M1,M2,M,PX1,PX2,PY1,PY2,M1X,M2X,M2Y,NPZ,AA,BB
	DOUBLE COMPLEX DLGM2X(2000),DLGM1X(2000),DLGM2Y(2000)
	DOUBLE COMPLEX DDM2X(2),DDM2Y(2),GU1,GU2,XX,YY,ZETA(2),S1,S2
	DOUBLE COMPLEX AN(20000),BN(20000),SU1(2000),SU2(2000),CCBACK
	DOUBLE PRECISION X,Y,NPY,AF,K1,K2,CHI(2000),PSI(2000)
	DOUBLE PRECISION CSCA,CEXT,NREAL,IX,IY,FD,FU(2003),MU(2000)
	DOUBLE PRECISION PI(2000,2000),TAU(2000,2000),Z,XY,XSUBI
	DOUBLE PRECISION D2BUCLE,SIGMA,BA,TH(2000),CBACK
	INTEGER NMAX,N,APROX,OPCION,K,T,P,PP,NMI
	COMMON/AK/X,Z,XY,D2BUCLE,M1,M2,SIGMA,P,PP,TH,OPCION,CONV                   
	EXTERNAL RBESSELR,DELOG
	IF (OPCION.EQ.1) THEN
		IX=X
		IY=IX+XSUBI
	ELSE IF (OPCION.EQ.2) THEN
		IX=XSUBI
		IY=IX+Z
	ELSE IF (OPCION.EQ.3) THEN
		IY=XSUBI                              
		IX=IY*D2BUCLE
	END IF
	IF(SIGMA.EQ.0) THEN 
		FD=1.0D0
	ELSE
		FD=(EXP(-0.5D0*(LOG(XSUBI/XY)/SIGMA)**2.0D0))/XSUBI
	END IF
        NMAX=INT(IY+4.0D0*IY**0.333333+2.0d0)
	M1X=M1*IX
	M2X=M2*IX
	M2Y=M2*IY
        M=M2/M1
	AF=1.0D-9
        IF(IX.NE.0.AND.M1.NE.M2) THEN
		CALL DELOG(M1X,NMAX,DLGM1X)
		IF(IX.NE.IY) THEN
			CALL DELOG(M2X,NMAX,DLGM2X)
			DDM2X(2)=-SIN(M2X)/COS(M2X)
			PX2=-DDM2X(2)
		END IF
	END IF
	CALL DELOG(M2Y,NMAX,DLGM2Y)
	CALL RBESSELR(IY,NMAX,CHI,PSI)
	DDM2Y(2)=-SIN(M2Y)/COS(M2Y)          
	PY2=-DDM2Y(2)
	ZETA(2)=(1.0D0,0.0D0)*CHI(1)-(0.0D0,1.0D0)*PSI(1)
	CSCA=0.0D0
	CEXT=0.0D0
	CBACK=0.0D0
	CCBACK=(0.0D0,0.0D0)
	IF (IX.EQ.0.OR.M1.EQ.M2) THEN
		APROX=1
	ELSE
		APROX=0
	END IF
	DO 1500 N=1,NMAX
	IF(APROX.EQ.1) THEN             
		XX=DLGM2Y(N+1)
		YY=DLGM2Y(N+1)
	ELSE IF (IX.EQ.IY) THEN
		XX=DLGM1X(N+1)*M
		YY=DLGM1X(N+1)/M
	ELSE
		NPZ=DFLOAT(N)/M2X
		DDM2X(1)=DDM2X(2)
		DDM2X(2)=1.0D0/(NPZ-DDM2X(1))-NPZ
		PX1=PX2
		PX2=PX1*(DDM2X(2)+NPZ)/(DLGM2X(N+1)+NPZ)
		NPZ=DFLOAT(N)/M2Y
		DDM2Y(1)=DDM2Y(2)
		DDM2Y(2)=1.0D0/(NPZ-DDM2Y(1))-NPZ
		PY1=PY2
		PY2=PY1*(DDM2Y(2)+NPZ)/(DLGM2Y(N+1)+NPZ)
		AA=(M*DLGM1X(N+1)-DLGM2X(N+1))/(M*DLGM1X(N+1)-DDM2X(2))
		AA=AA*PX2
		BB=(DLGM1X(N+1)-M*DLGM2X(N+1))/(DLGM1X(N+1)-M*DDM2X(2))
		BB=BB*PX2
		K1=MAX(ABS(AA),ABS(BB))
		K2=AF*ABS(PY2)*MIN(ABS(DLGM2Y(N+1)/DDM2Y(2)),1.0D0)
                IF (K1.LT.K2) THEN
                        XX=DLGM2Y(N+1)
                        YY=DLGM2Y(N+1)
                        APROX=1
                ELSE
			GU1=DLGM2Y(N+1)*PY2
			XX=(GU1-DDM2Y(2)*AA)/(PY2-AA)
			YY=(GU1-DDM2Y(2)*BB)/(PY2-BB)
                END IF
		CONTINUE
	END IF
		NPY=DFLOAT(N)/IY
		GU1=XX/M2+NPY
		GU2=YY*M2+NPY
		ZETA(1)=ZETA(2)
		ZETA(2)=(1.0D0,0.0D0)*CHI(N+1)-(0.0D0,1.0D0)*PSI(N+1)
		AN(N)=(GU1*CHI(N+1)-CHI(N))/(GU1*ZETA(2)-ZETA(1))
		BN(N)=(GU2*CHI(N+1)-CHI(N))/(GU2*ZETA(2)-ZETA(1))
		NREAL=DFLOAT(N+N)+1.0D0
		CEXT=CEXT+NREAL*DREAL(AN(N)+BN(N))
		CSCA=CSCA+NREAL*(ABS(AN(N))**2+ABS(BN(N))**2)
		CCBACK=CCBACK+NREAL*(AN(N)-BN(N))*(-1.0D0)**N
1500    CONTINUE
	CBACK=(ABS(CCBACK)**2)
	IF(P.NE.0) THEN
	DO K=1,P
	S1=(0.0D0,0.0D0)
	S2=(0.0D0,0.0D0)
	DO T=1,NMAX
		IF(PI(K,1).NE.1) THEN
			MU(K)=COS(TH(K))
			PI(K,1)=1.0D0
			PI(K,2)=3.0D0*MU(K)
			TAU(K,1)=MU(K)
			TAU(K,2)=2.0D0*MU(K)*PI(K,2)-3.0D0
			NMI=2
		END IF
		RT=DFLOAT(T)                
		IF(T.GT.NMI) THEN
	PI(K,T)=((2.0D0*RT-1.0D0)*MU(K)*PI(K,T-1)-RT*PI(K,T-2))/(RT-1.0D0)
	TAU(K,T)=RT*MU(K)*PI(K,T)-(RT+1.0D0)*PI(K,T-1)
		END IF
		BA=(2.0D0*RT+1.0D0)/(RT*(RT+1.0D0))
		S1=S1+(AN(T)*PI(K,T)+BN(T)*TAU(K,T))*BA
          S2=S2+(BN(T)*PI(K,T)+AN(T)*TAU(K,T))*BA                             
	END DO
700     CONTINUE
	FU(K)=(ABS(S1)*ABS(S1)+abs(s2)*abs(s2))*0.5d0*FD
	SU1(K)=S1
	SU2(K)=S2
	END DO
	IF(NMAX.GT.NMI) NMI=NMAX
	END IF
	FU(P+1)=CEXT*FD
	FU(P+2)=CSCA*FD
	FU(P+3)=CBACK*FD
	RETURN
	END
C
C
C
C
	SUBROUTINE RBESSELR(Y,NMAX,CHI,PSI)
	DOUBLE PRECISION Y,CHI(20000),PSI(20000),JJ(20200),JZ,QA
	INTEGER NMAX,N1,TT
	IF(NMAX.GT.150) THEN
		N1=IDNINT(1.1D0*NMAX)
	ELSE
		N1=NMAX+15
	END IF
	JJ(N1)=0.0D0
	JJ(N1-1)=1.0D-35
	DO TT=N1-2,1,-1
		JJ(TT)=(DFLOAT(TT+TT)+3.0D0)*JJ(TT+1)/Y-JJ(TT+2)
	END DO
	JZ=3.0D0*JJ(1)/Y-JJ(2)
	QA=(SIN(Y))/JZ
	CHI(1)=SIN(Y)
	CHI(2)=JJ(1)*QA
	PSI(1)=COS(Y)
	PSI(2)=PSI(1)/Y+SIN(Y)
	DO TT=2,NMAX+1
		CHI(TT+1)=JJ(TT)*QA
		PSI(TT+1)=(DFLOAT(TT+TT)-1.0D0)*PSI(TT)/Y-PSI(TT-1)
	END DO
	CONTINUE
	RETURN
	END
C
C
C
C
	SUBROUTINE DELOG(Z,NMAX,DN)
	DOUBLE COMPLEX Z,TK,DN(20000)
	INTEGER NMAX,N1,TT
        N1=MAX(INT(ABS(Z)),NMAX)
	IF(N1.GT.150) THEN
		N1=IDNINT(1.1D0*N1)
	ELSE
		N1=N1+15
	END IF
	DN(N1)=(0.0D0,0.0D0)
	DO TT=N1-1,1,-1
		TK=DFLOAT(TT)/Z
		DN(TT)=TK-1.0D0/(DN(TT+1)+TK)
	END DO
	CONTINUE
	RETURN
	END