Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

cfffv11.F @ 201

Last change on this file since 201 was 175, checked in by slebonnois, 14 years ago

S.LEBONNOIS:

Revision majeure de la physique Titan => ajout des nuages version 10 bins (Jeremie Burgalat) Cette version reste a tester mais avec clouds=0, on reste sur l'ancienne.
Quelques ajouts dans la doc.

File size: 36.3 KB

Rev	Line
[175]	1	subroutine cfffv11(lambda,xn,xk,rad,DFS,NB,fsca,fext,fabs,fg)
	2	*
	3	* NEW VERSION MARCH, 31, 1999.
	4	* WITHOUT MAXWELL-GARNETT APPROACH FOR IM{M} > 0.1
	5	*
	6
	7
	8	parameter (nang=91)
	9
	10	complex mtest,refrel,s1(2000),s2(2000)
	11	real rad,lambda,s11(2000),theta(10000)
	12	real s110u(2000),s111u(2000)
	13	real s110d(2000),s111d(2000)
	14	real pol(2000),pp(2,2000),strS(2000)
	15	real NB
	16	real xxn
	17
	18	* COMMON WITH INTMIE
	19
	20	common/angle11/theta
	21
	22	* COMMON WITH CORRINT11
	23
	24	common/pack1/cjup1,cju,cjup1_,cju_,xrap,xechj,xechjp1
	25	common/pack2/cjdp1,cjd,cjdp1_,cjd_,xku,xkd,jindex
	26	common/pack3/xexp5,xexp6
	27
	28	* COMMON WITH FFFV11
	29
	30	DATA IND /0/
	31
	32
	33	if(nang*2.gt.2000)
	34	& stop 'INCREASE THE SIZE OF THE ARRAYS IN CFFFV11'
	35
	36	c refrel=(xn,xk)
	37	refrel=cmplx(xn,xk)
	38	nindex=int(xn*100.)
	39	x=2.3.14159265rad/lambda
	40	dang=1.570796327/dfloat(nang-1)
	41	pi=3.14159265
	42	NC=20
	43	DF=1.99999
	44	alpha=1.1
	45	alphaS=(1.4-1.1)/(2.5-2.)*DFS-.1 !APPROXIMATED!
	46	!EXACT FOR D=2 AND 2.5!
	47
	48	******************************************************************
	49	* STEP # 1
	50	******************************************************************
	51	* CALL THE ROUTINE THAT GIVE THE SCALING FACTOR'S PSI_s AND PSI_e
	52	* FOR EACH BOUND OF THE INTERVAL nX(j) < n*X < nX(j+1) AND
	53	* xku < xk < xkd.
	54	* THESE FACTORS ARE GIVEN ASSUMING : N=20, DF=2 AND KNOWING n
	55	*
	56	* NB: HERE PSI ARE NOTED cju,cjd,cjup1,cjdp1.....
	57	******************************************************************
	58
	59	call corrint11(x,xn,xk,NC)
	60
	61
	62	******************************************************************
	63	* STEP # 2
	64	******************************************************************
	65	* FOR EACH OF THE FOUR "POINTS" (nX,k), COMPUTE THE MIE CROSS-SECTIONS
	66	* AND PHASE FUNCTION FOR A SPHERE AS ONE MONOMER.
	67	* DERIVE, THANKS TO THE PSI, THE MONOMER CROSS-SECTIONS
	68	* COMPUTE ALSO THE ACTUAL MIE CROSS SECTION (I.E FOR THE ACTUAL
	69	* VALUES OF n,X and k)
	70	******************************************************************
	71
	72	* STEP #2.1
	73	******************************************
	74	* if Mie & monomer in Rayleigh scattering:
	75	******************************************
	76	* the phase function shape does not depends
	77	* on x neither the PSI's: for each values of
	78	* k, only one computation is to be done
	79	*****************************************
	80
	81	IF (x*xn.le.0.085) THEN
	82
	83
	84	* pour k=xku
	85	* # x=xech(j)
	86	* # x=xech(j+1)
	87
	88	c refrel=(xn,xku)
	89	refrel=cmplx(xn,xku)
	90
	91	call intmie(x,refrel,nang,s1,s2,qext_,qsca_)
	92	qsca1u=alog10(qsca_*400.)
	93	qext1u=alog10(qext_*20.)
	94	qsca2u=alog10(qsca_*400.)
	95	qext2u=alog10(qext_*20.)
	96
	97	* pour k=xkd
	98	* # x=xech(j)
	99	* # x=xech(j+1)
	100
	101	c refrel=(xn,xkd)
	102	refrel=cmplx(xn,xkd)
	103
	104	call intmie(x,refrel,nang,s1,s2,qext_,qsca_)
	105	qsca1d=alog10(qsca_*400.)
	106	qext1d=alog10(qext_*20.)
	107	qsca2d=alog10(qsca_*400.)
	108	qext2d=alog10(qext_*20.)
	109
	110	xrap=0.
	111
	112	* Mie exact
	113	*
	114
	115	c refrel=(xn,xk)
	116	refrel=cmplx(xn,xk)
	117
	118	call intmie(x,refrel,nang,s1,s2,qext_,qsca_)
	119	do j=1,2*nang-1
	120	s11(j)=cabs(s2(j))cabs(s2(j))+cabs(s1(j))cabs(s1(j))
	121	s11(j)=s11(j)/(2.pix*2.qsca_)
	122	enddo
	123
	124
	125	ELSE
	126
	127	* STEP #2.2
	128	************************************************
	129	* if Mie & monomer in the range X about 1 to 5 *
	130	************************************************
	131
	132	* pour k=xkd
	133	* # x=xech(j)
	134	* # x=xech(j+1)
	135
	136
	137	qsca1d=alog10(cjd)
	138	qext1d=alog10(cjd_)
	139
	140	qsca2d=alog10(cjdp1)
	141	qext2d=alog10(cjdp1_)
	142
	143	* pour k=xku
	144	* # x=xech(j)
	145	* # x=xech(j+1)
	146
	147
	148	qsca1u=alog10(cju)
	149	qext1u=alog10(cju_)
	150
	151	qsca2u=alog10(cjup1)
	152	qext2u=alog10(cjup1_)
	153
	154
	155	* Mie exact:
	156	*
	157
	158	c refrel=(xn,xk)
	159	refrel=cmplx(xn,xk)
	160	call intmie(x,refrel,nang,s1,s2,qext_,qsca_)
	161	do j=1,2*nang-1
	162	s11(j)=cabs(s2(j))cabs(s2(j))+cabs(s1(j))cabs(s1(j))
	163	enddo
	164
	165	ENDIF
	166
	167	* THIS IS MIE:
	168
	169	qabs_ =qext_ -qsca_
	170
	171	* THIS IS FRACTAL 20 MONOMERS
	172	* d
	173	qabs2d=alog10(10.qext2d-10.qsca2d)
	174	qabs1d=alog10(10.qext1d-10.qsca1d)
	175	* u
	176	qabs2u=alog10(10.qext2u-10.qsca2u)
	177	qabs1u=alog10(10.qext1u-10.qsca1u)
	178
	179
	180	***********************************************************************
	181	* QUICK FIX FOR NEGATIVE ABSORPTION WHEN APPEARS [OUT OF n RANGE]
	182	**********************************************************************
	183	GOTO 505
	184	* This avoid a "crash" if one among four of the qsca is larger than qext
	185	* Averaging the w=qsca/qext, and applied to the "wrong" qsca
	186	* NOTE this occurs only if real refractive index > 2
	187
	188	w1=10.qsca2d/10.qext2d
	189	w2=10.qsca1d/10.qext1d
	190	w3=10.qsca2u/10.qext2u
	191	w4=10.qsca1u/10.qext1u
	192
	193
	194	IF(10.qext2d.lt.10.qsca2d)
	195	& qsca2d=alog10((10.*qext2d)(w2+w3+w4)/3.)
	196	IF(10.qext1d.lt.10.qsca1d)
	197	& qsca1d=alog10((10.*qext1d)(w1+w3+w4)/3.)
	198	IF(10.qext2u.lt.10.qsca2u)
	199	& qsca2u=alog10((10.*qext2u)(w2+w1+w4)/3.)
	200	IF(10.qext1u.lt.10.qsca1u)
	201	& qsca1u=alog10((10.*qext1u)(w2+w3+w1)/3.)
	202
	203	qabs2d=alog10(10.qext2d-10.qsca2d)
	204	qabs1d=alog10(10.qext1d-10.qsca1d)
	205	qabs2u=alog10(10.qext2u-10.qsca2u)
	206	qabs1u=alog10(10.qext1u-10.qsca1u)
	207	505 continue
	208
	209	************************************************
	210	* STEP # 5
	211	************************************************
	212	* INTERPOLATION OVER THE TWO VARIABLES n*X AND K
	213	************************************************
	214
	215	* X*n VARIABLE/ interpolation log-log.
	216
	217	sextu=((qext2u-qext1u)*xrap+qext1u)
	218	sscau=((qsca2u-qsca1u)*xrap+qsca1u)
	219	sabsu=((qabs2u-qabs1u)*xrap+qabs1u)
	220	sextd=((qext2d-qext1d)*xrap+qext1d)
	221	sscad=((qsca2d-qsca1d)*xrap+qsca1d)
	222	sabsd=((qabs2d-qabs1d)*xrap+qabs1d)
	223
	224	* K VARIABLE/ interpolation log-log
	225	xki=xk
	226	if (xk.lt.1.e-2) xki=1.e-2
	227	if (xk.gt.1.) xki=1.
	228	rki=-alog10(xki)
	229	rku=-alog10(xku)
	230	rkd=-alog10(xkd)
	231	ssca=(sscau-sscad)*(rki-rkd)+sscad
	232	sext=(sextu-sextd)*(rki-rkd)+sextd
	233	sabs=(sabsu-sabsd)*(rki-rkd)+sabsd
	234
	235	sscap=(sscau-sscad)*(rki-rkd)+sscad
	236	sextp=(sextu-sextd)*(rki-rkd)+sextd
	237	sabsp=(sabsu-sabsd)*(rki-rkd)+sabsd
	238
	239	* storage Mie down and up
	240	c refrel=(xn,1.)
	241	refrel=cmplx(xn,1.)
	242	call intmie(x,refrel,nang,s1,s2,qextmd_,qscamd_)
	243	qabsmd_=qextmd_-qscamd_
	244
	245	c refrel=(xn,.1)
	246	refrel=cmplx(xn,.1)
	247	call intmie(x,refrel,nang,s1,s2,qextmu_,qscamu_)
	248	qabsmu_=qextmu_-qscamu_
	249
	250	******************************************************
	251	* STEP # 5.1 :
	252	******************************************************
	253	* SPECIAL CASE FOR LARGE k (>.1) : MAXWELL GARNETT
	254	******************************************************
	255
	256
	257	if (xk.gt.0.1) then ! THIS CONCERNS THE LARGE IMAGINARY REFRACTIVE INDEXES
	258
	259	*Prepare Maxwell Garnett approx. (see B&H)
	260	*-----------------------------------------
	261
	262	xbulk=x(NC1.)**.3333333
	263	xsphe=x(NC1.)**.5
	264	frac=(xbulk/xsphe)**3.
	265	xkff=xk*frac
	266	xnff=1.+(xn-1.)*frac
	267
	268	* For small x, aggregate scattering= N**2 monmer scatt.:
	269	* -----------------------------------------------------
	270
	271	mtest=cmplx(xn,xk)
	272	call intmie(x,mtest,nang,s1,s2,qe,qs)
	273	sscaR=alog10(qsNC(xsphe/x)**2.)
	274
	275
	276	* For small x, aggregate abs. is proportionnal to N
	277	* -----------------------------------------------------
	278	mtest=cmplx(xn,xk)
	279	call intmie(x,mtest,nang,s1,s2,qe,qs)
	280	sabsR=alog10((qe-qs)*20.)
	281
	282
	283	endif
	284
	285
	286
	287	*****************************************************************
	288	* FINALLY, WE GET THE VALUES FOR AN AGGREGATE OF 20 MONOMERS
	289	* FOR THE ACTUAL X,k, AND n
	290	*****************************************************************
	291
	292
	293	ssca=10.**ssca
	294	sext=10.**sext
	295	sabs=10.**sabs
	296
	297	* Sharp cross over between Rayleigh
	298	xup=0.3 ! \|<----------
	299	xdo=0.3 ! --------->\|
	300
	301	if (xk.gt.0.1) then ! large imaginary refractive index
	302
	303	if (x*xn.ge.xup) then ! large size parameter
	304	sext=sabs+ssca
	305	scal1=(10.**sabsd)/qabsmd_
	306	scal2=(10.**sabsu)/qabsmu_
	307	scalx=(alog10(xk)+1.)*(scal1-scal2)+scal2
	308	sabs=qabs_*scalx
	309	ssca=sext-sabs
	310	endif
	311
	312	if (x*xn.le.xdo) then ! PURE RAYLEIGH
	313	sabs=10.**sabsR
	314	ssca=10.**sscaR
	315	sext=sabs+ssca
	316	endif
	317
	318
	319	else
	320
	321	sext=sabs+ssca !self-consistent ext sca abs set
	322
	323	endif
	324
	325
	326	*****************************************************************
	327	* STEP # 6
	328	*****************************************************************
	329	* NOW, USE USE THE N-LAWS TO COMPUTE THE Q(N) FROM Q(N=20) AND
	330	* THE Q(MIE) [EQUATIONS (9) AND (10)].
	331	*****************************************************************
	332
	333	g0=0.
	334	g1=0.
	335	g2=0.
	336
	337	surf=rad*2.3.1415926*1.e18
	338
	339
	340	* STEP # 6.1
	341	*****************************************************************
	342	* FIRST; COMPUTE THE STRUCTURE TERM THAT IS ONLY USED TO COMPUTE
	343	* THE ASYMETRY PARAMETER THAT DEPENDS ON THE SHAPE OF THE PHASE
	344	* FUNCTION.
	345	*****************************************************************
	346
	347
	348	do 367 j=1,2*nang-1
	349	z=sin(theta(j)/2.)
	350	xx_=alphaSx(NB1.)*(1./DFS)
	351	call structure(NB,xx_,z,str,DFS)
	352	if (z.eq.0.) str=(NB1.)*2.
	353	strS(j)=str
	354	367 continue
	355
	356
	357	do 365 j=1,2*nang-2,2
	358
	359	g1=2dang/6(s11(j)sin(theta(j))strS(j)
	360	& +s11(j+2)sin(theta(j+2))strS(j+2)
	361	& +4s11(j+1)sin(theta(j+1))*strS(j+1))
	362	& 23.141592+g1
	363
	364	g2=2dang/6(s11(j)sin(theta(j))cos(theta(j))*strS(j)
	365	& +s11(j+2)sin(theta(j+2))cos(theta(j+2))*strS(j+2)
	366	& +4s11(j+1)sin(theta(j+1))cos(theta(j+1))strS(j+1))
	367	& 2.3.141592+g2
	368
	369	365 continue
	370
	371
	372	* STEP # 6.2
	373	*****************************************************************
	374	* USE THE EQUATIONS (10) AND (9) FOR MEDIUM RANGE OR THE RAYLEIGH
	375	* RANGE.
	376	*****************************************************************
	377
	378
	379	* 1: FOR MEDIUM X RANGE: 1< X < 10 to 50 [EQ 10]
	380	* 2: FOR VERY LARGE X RANGE: X >> 50 [EQ 11]
	381	* N^a LOG(N^b) IS REPLACED BY A PURE N^a LAW. THIS OCCURS
	382	* WHEN N>50 FOR ABSORPTION AND WHEN N>25/x FOR SCATTERING,
	383	* AND CORRESPOND TO THE GEOMETRIC RANGE. THE EXPONENT a IS
	384	* SET BY EMPIRIC CONSIDERATIONS (SEE PAPER).
	385
	386	* XEXP5 AND XEXP6 ARE THE EXPONENT FACTOR USED IN EQ 10
	387	* XEXP1 AND XEXP2 ARE THE EXPONENT FACTOR USED VERY LARGE BEHAVIOUR
	388
	389	xkf=1.
	390	if (xk.le.1.e-2) xkf=xk/1.e-2
	391
	392
	393	rho=(NC1.)*xexp5
	394	rho_=(NB1.)*xexp5
	395
	396	* EQ 10 RESULT FOR ABSORPTION:
	397	*-----------------------------
	398
	399	sabs1=(sabsxkf-qabs_)/alog10(20.)alog10(NB*1.)
	400	& *rho_/rho+qabs_
	401
	402
	403	* NB: IF NB>NTA OR NB>NTS , EFFICIENCY FACTOR FOLLOW A POWER LAW
	404	* BEHAVIOUR. ELSE, EQUATION 10 IS VALIDE.
	405
	406	NTA=50
	407
	408	if(NB.gt.NTA) then
	409	rhox=(NTA1.)*xexp5
	410
	411	* COMPUTE sabs1 FOR NTA TO START THE POWER LAW:
	412	sabs1=(sabsxkf-qabs_)/alog10(20.)alog10(NTA*1.)
	413	& *rhox/rho+qabs_
	414
	415	* COMPUTE xexp1 EXPONENT FACTOR:
	416	xexp1=.85
	417	if(xxn.lt.2.2) xexp1=(.85-.96)1.42(xxn-1.5)+.96
	418	if(x*xn.lt.1.5) xexp1=.96
	419
	420	* COMPUTE sabs1 VALUE FOR LARGE X:
	421	sabs1=sabs1(NB1./(NTA1.))*xexp1
	422	endif
	423
	424	* HERE SABS1 MEDIUM RANGE IS KNOWN
	425
	426	rho=(NC1.)*xexp6
	427	rho_=(NB1.)*xexp6
	428
	429
	430	* EQ 10 RESULT FOR SCATTERING:
	431	*-----------------------------
	432
	433	ssca1=(ssca-qsca_)/alog10(20.)alog10(NB1.)
	434	& *rho_/rho+qsca_
	435
	436	NTS=int((20./x)**(1./.65))
	437
	438	* NB: IF NB>NTA OR NB>NTS , EFFICIENCY FACTOR FOLLOW A POWER LAW
	439	* BEHAVIOUR. ELSE, EQUATION 10 IS VALIDE.
	440
	441	if(NB.gt.NTS) then
	442	rhox=(NTS1.)*xexp6
	443
	444	* COMPUTE sabs1 FOR NTS TO START THE POWER LAW:
	445	ssca1=(ssca-qsca_)/alog10(20.)alog10(NTS1.)
	446	& *rhox/rho+qsca_
	447
	448	* COMPUTE xexp2 EXPONENT FACTOR:
	449	xexp2=.95
	450	if(x.lt.1.0) xexp2=(.95-1.1)2(x-.5)+1.1
	451
	452	* COMPUTE sabs1 VALUE FOR LARGE X:
	453	ssca1=ssca1(NB1./(NTS1.))*xexp2
	454	endif
	455
	456	* HERE SSCA1 MEDIUM RANGE IS KNOWN
	457
	458
	459	sext1=ssca1+sabs1
	460
	461
	462
	463
	464
	465	* 3: FOR THE RAYLEIGH RANGE [EQ 9]
	466
	467	ssca2=(alog10(ssca)-alog10(qsca_))/alog10(20.)
	468	& alog10(NB1.)+alog10(qsca_)
	469	sext2=(alog10(sext)-alog10(qext_))/alog10(20.)
	470	& alog10(NB1.)+alog10(qext_)
	471	sabs2=(alog10(sabs*xkf)-alog10(qabs_))/alog10(20.)
	472	& alog10(NB1.)+alog10(qabs_)
	473
	474
	475
	476	************************************************
	477	* STEP # 7
	478	************************************************
	479	* THE CROSS OVER BETWEEN RAYLEYGH AND MEDIUM RANGE
	480	************************************************
	481
	482	rho_=(NB1.)*.66666
	483	* 1< X < 50.
	484	sext1=alog10(sext1/rho_)
	485	ssca1=alog10(ssca1/rho_)
	486	sabs1=alog10(sabs1/rho_)
	487	* RAYLEIGH
	488	sext2=alog10((10.**sext2)/rho_)
	489	ssca2=alog10((10.**ssca2)/rho_)
	490	sabs2=alog10((10.**sabs2)/rho_)
	491
	492
	493
	494
	495
	496
	497	* 7.1: SCATTERING AND EXTINCTION + CROSS OVER
	498	************************************************
	499
	500
	501	* RAYLEIGH RANGE:
	502	fsca=10.**(ssca2)
	503	fext=10.**(sext2)
	504	fg=g2/g1
	505
	506	xx=.5alphax(NB1.)**.33333
	507	xup=10.
	508	xdo=.1
	509
	510	f=alog10(xx/xdo)/alog10(xup/xdo)
	511	g=1.-f
	512
	513
	514	* MEDIUM X
	515	if (xx.le.xup.and.xx.ge.xdo) then
	516	fsca=10.*(fssca1+g*ssca2)
	517	fext=10.*(fsext1+g*sext2)
	518	fg=g2/g1
	519	endif
	520
	521	* LARGE X
	522	if (xx.gt.xup) then
	523	fsca=10.**(ssca1)
	524	fext=10.**(sext1)
	525	fg=g2/g1
	526	endif
	527
	528	* 7.2: ABSORPTION + CROSS OVER
	529	************************************************
	530
	531	* RAYLEIGH RANGE
	532
	533	fabs=10.**(sabs2)
	534
	535	xx=.5alphax(NB1.)**.5
	536	f=alog10(xx/xdo)/alog10(xup/xdo)
	537	g=1.-f
	538
	539	* MEDIUM X
	540	if (xx.le.xup.and.xx.ge.xdo) then
	541	fabs=10.*(fsabs1+g*sabs2)
	542	endif
	543
	544	* LARGE X
	545	if (xx.gt.xup) then
	546	fabs=10.**(sabs1)
	547	endif
	548
	549	fext=fsca+fabs
	550
	551
	552	***********************************************
	553	* STEP # 8
	554	***********************************************
	555	* DISPLAY THE RESULTS : PHASE FUNCTION OR
	556	* EFFICIENCY COEFFICIENTS
	557	***********************************************
	558
	559	IF (IND.eq.1) THEN
	560
	561	* 8.1: CROSS SECTIONS & PHASE FUNCTION
	562	***********************************************
	563
	564	tot =0.
	565	scoef=NB*(2./3.)rad*2.3.1415926
	566
	567	do 366 j=1,2*nang-1
	568
	569	pp(1,j)=theta(j)*180./3.1415926
	570	pp(2,j)=s11(j)23.1415936strS(j)1./g1
	571	366 continue
	572	*
	573	ELSE
	574
	575	* 8.2: EFFICIENCY COEFFICIENTS
	576	***********************************************
	577
	578	scoef=NB*(2./3.)rad*2.3.1415926 ! METERS^2
	579	fsca=fsca*scoef
	580	fext=fext*scoef
	581	fabs=fabs*scoef
	582
	583	c new fashion...
	584
	585	fabs=(qext_-qsca_)NB/(NB1.)*(2./3.)scoef
	586	fext=fsca+fabs
	587
	588
	589	ENDIF
	590
	591	return
	592	end
	593
	594	c------------------------------------------------------------
	595	subroutine intmie(x,refrel,nang,s1,s2,
	596	& qext,qsca)
	597	**********************************************************
	598	* THIS ROUTINE COMES FROM BORHEN AND HUFFMAN BOOK:
	599	* "ABSORPTION AND SCATTERING OF LIGHT BY SMALL PARTICLES"
	600	* WILEY INTERSCIENCE PUBLICATION, 1983
	601	**********************************************************
	602	common/angle11/theta
	603
	604	real amu(10000),theta(10000),pi(10000)
	605	real tau(10000),pi0(10000),pi1(10000)
	606	complex d(300000),y,refrel,xi,xi0,xi1,an,bn,s1(20000),s2(20000)
	607	complex s_(20000)
	608	double precision psi0,psi1,psi,dn,dx
	609	dx=x !dx en double precision ....
	610	y=x*refrel
	611
	612	xstop=x+4x*.3333+2.
	613	nstop=xstop
	614	ymod=cabs(y)
	615	nmx=amax1(xstop,ymod)+15
	616	c print*,nmx,xstop,nstop,ymod
	617	dang=1.570796327/dfloat(nang-1)
	618	do 555 j=1,2*nang-1
	619	theta(j)=(dfloat(j)-1.)*dang
	620	555 amu(j)=cos(theta(j))
	621	c d(nmx)=(0.,0.)
	622	d(nmx)=cmplx(0.,0.)
	623	nn=nmx-1
	624
	625	do 120 n=1,nn
	626	rn=nmx-n+1
	627	d(nmx-n)=(rn/y)-(1./(d(nmx-n+1)+rn/y))
	628	120 continue
	629
	630	do 666 j=1,nang
	631	pi0(j)=0. ! fonction de legendre
	632	666 pi1(j)=1.
	633
	634	nn=2*nang-1
	635
	636	do 777 j=1,nn
	637	c s_(j)=(0.,0.)
	638	c s1(j)=(0.,0.)
	639	c777 s2(j)=(0.,0.)
	640	s_(j)=cmplx(0.,0.)
	641	s1(j)=cmplx(0.,0.)
	642	777 s2(j)=cmplx(0.,0.)
	643	psi0=dcos(dx) !initialisation des fonctions de Bessel
	644	psi1=dsin(dx)
	645	chi0=-sin(x)
	646	chi1=cos(x)
	647	apsi0=psi0 !psi en double prec. et apsi en simple
	648	apsi1=psi1
	649	c xi0=(apsi0,-chi0)
	650	c xi1=(apsi1,-chi1)
	651	xi0=cmplx(apsi0,-chi0)
	652	xi1=cmplx(apsi1,-chi1)
	653	qsca=0.
	654	qsca_=0.
	655	n=1
	656	c ***********debut de l'iteration sur n ***********
	657	200 dn=n
	658	rn=n
	659	fn=(2.rn+1.)/(rn(rn+1.))
	660
	661	psi=(2.dn-1.)psi1/dx-psi0 ! calcul des fct de Bessel
	662	chi=(2.rn-1.)chi1/x-chi0 ! relation recurrente a 2 niveaux
	663	apsi=psi
	664	c xi=(apsi,-chi)
	665	xi=cmplx(apsi,-chi)
	666	an=(d(n)/refrel+rn/x)*apsi-apsi1
	667	an=an/((d(n)/refrel+rn/x)*xi-xi1)
	668	bn=(refreld(n)+rn/x)apsi-apsi1
	669	bn=bn/((refreld(n)+rn/x)xi-xi1)
	670	qsca=qsca+(2rn+1.)(cabs(an)cabs(an)+cabs(bn)cabs(bn))
	671
	672	c *************debut de la boucle sur les angles*****
	673	do 789 j=1,nang
	674	jj=2*nang-j
	675	pi(j)=pi1(j) !
	676	tau(j)=rnamu(j)pi(j)-(rn+1.)*pi0(j) ! fonction de legendre
	677	s1(j)=s1(j)+fn(anpi(j)+bn*tau(j))
	678	s2(j)=s2(j)+fn(antau(j)+bn*pi(j))
	679	p=(-1)**(n-1)
	680	t=(-1)**n
	681	if (j.eq.jj) goto 789
	682	s1(jj)=s1(jj)+fn(anpi(j)p+bntau(j)*t)
	683	s2(jj)=s2(jj)+fn(antau(j)t+bnpi(j)*p)
	684	789 continue
	685	psi0=psi1
	686	psi1=psi
	687	apsi1=psi1 ! double prec=simple
	688	chi0=chi1
	689	chi1=chi
	690	c xi1=(apsi1,-chi1)
	691	xi1=cmplx(apsi1,-chi1)
	692	n=n+1
	693	rn=n
	694	do 999 j=1,nang
	695	pi1(j)=((2.rn-1.)/(rn-1.))amu(j)*pi(j)
	696	pi1(j)=pi1(j)-rn*pi0(j)/(rn-1.)
	697	999 pi0(j)=pi(j)
	698	if (n-1-nstop) 200,300,300
	699	300 qsca=(2./(xx))qsca
	700	qext=(4./(xx))real(s1(1))
	701	qabs=qext-qsca
	702
	703
	704	return
	705	end
	706
	707
	708	subroutine corrint11(x,xn,xk,NB)
	709	parameter (nech=28)
	710	parameter (nex=6)
	711	real xech(nech)
	712	real A0(2nech), B0(2nech), C0(2nech), D0(2nech)
	713	real A1(2nech), B1(2nech), C1(2nech), D1(2nech)
	714	real A2(2nech), B2(2nech), C2(2nech), D2(2nech)
	715	real A3(2nech), B3(2nech), C3(2nech), D3(2nech)
	716	real x,xn,xk,correct,correct_
	717	integer NB,ifirst
	718	real ww1(nex),ww2(nex)
	719	real xx1(nex),xx2(nex)
	720	real yy1(nex),yy2(nex)
	721	real zz1(nex),zz2(nex)
	722	real xldo(nex)
	723	real xexp5,xexp6
	724	common/pack1/cjup1_,cju_,cjup1,cju,xrapl,xechj,xechjp1
	725	common/pack2/cjdp1_,cjd_,cjdp1,cjd,xku,xkd,jindex
	726	common/pack3/xexp5,xexp6
	727	* * ATTENTION: ORDRE INVERSE DANS cfffpf, cfffcs,... *
	728	DATA xech/0.05,0.1,0.25,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,
	729	& 1.4,1.5,1.75,2.0,2.25,2.5,2.75,3.0,3.25,3.5,3.75,4.0,4.25,
	730	& 4.5,4.75,5.0/
	731	DATA A0/ .1117E+00, .2284E+00, .4473E+00,-.3351E+00,-.7323E+00,
	732	& -.8070E+00,-.5968E+00,-.2956E+00, .1171E-01, .2515E+00,
	733	& .3607E+00, .3653E+00, .2658E+00, .1145E+00,-.1505E+00,
	734	& .1982E-01, .8496E-01,-.5869E-01,-.7623E-01,-.5410E-02,
	735	& -.2850E-01,-.5501E-01,-.5267E-01,-.2317E-01,-.3344E-01,
	736	& -.3936E-01,-.3094E-01,-.3099E-01,
	737	& .6136E-03, .7863E-02, .9823E-01, .6704E-01,-.8727E-01,
	738	& -.1741E+00,-.1402E+00,-.3270E-01, .1074E+00, .2397E+00,
	739	& .3254E+00, .3654E+00, .3495E+00, .2869E+00, .7161E-01,
	740	& .6914E-01, .1319E+00, .7578E-01, .2340E-01, .4031E-01,
	741	& .4202E-01, .2560E-01, .1679E-01, .2965E-01, .3108E-01,
	742	& .2983E-01, .3023E-01, .2988E-01/
	743	DATA B0/-.5034E+00,-.1033E+01,-.2186E+01, .3670E+00, .2014E+01,
	744	& .2646E+01, .2279E+01, .1496E+01, .5713E+00,-.2137E+00,
	745	& -.6203E+00,-.7014E+00,-.4274E+00, .3641E-01, .8890E+00,
	746	& .3026E+00, .4121E-01, .4774E+00, .5016E+00, .2491E+00,
	747	& .3097E+00, .3754E+00, .3536E+00, .2500E+00, .2734E+00,
	748	& .2838E+00, .2453E+00, .2429E+00,
	749	& -.2493E-02,-.3224E-01,-.4251E+00,-.5115E+00, .7309E-01,
	750	& .5198E+00, .5940E+00, .3982E+00, .3812E-01,-.3514E+00,
	751	& -.6344E+00,-.7949E+00,-.7802E+00,-.6081E+00, .5837E-01,
	752	& .4640E-01,-.1963E+00,-.4932E-01, .9157E-01, .1454E-01,
	753	& -.6767E-02, .3216E-01, .4963E-01,-.9453E-04,-.1184E-01,
	754	& -.1405E-01,-.2189E-01,-.2395E-01/
	755	DATA C0/ .6106E+00, .1257E+01, .2904E+01, .1547E+01, .1531E+00,
	756	& -.5216E+00,-.4001E+00, .1061E+00, .7939E+00, .1421E+01,
	757	& .1776E+01, .1883E+01, .1693E+01, .1334E+01, .6370E+00,
	758	& .1107E+01, .1329E+01, .9812E+00, .9691E+00, .1171E+01,
	759	& .1119E+01, .1070E+01, .1086E+01, .1164E+01, .1145E+01,
	760	& .1135E+01, .1167E+01, .1164E+01,
	761	& .2614E-02, .3429E-01, .4919E+00, .1051E+01, .6861E+00,
	762	& .3474E+00, .2567E+00, .3720E+00, .6355E+00, .9464E+00,
	763	& .1195E+01, .1361E+01, .1388E+01, .1283E+01, .7888E+00,
	764	& .8176E+00, .1038E+01, .9433E+00, .8490E+00, .9234E+00,
	765	& .9487E+00, .9256E+00, .9171E+00, .9601E+00, .9731E+00,
	766	& .9774E+00, .9868E+00, .9892E+00/
	767	DATA A1/ .8917E-02, .1683E-01, .1060E-01,-.2180E+00,-.4334E+00,
	768	& -.6977E+00,-.9510E+00,-.1114E+01,-.1154E+01,-.1186E+01,
	769	& -.1170E+01,-.1184E+01,-.1214E+01,-.1243E+01,-.1324E+01,
	770	& -.1570E+01,-.1304E+01,-.7215E+00,-.9031E+00,-.1164E+01,
	771	& -.5731E+00,-.4235E+00,-.8438E+00,-.3892E+00, .3259E+00,
	772	& .2680E+00, .6864E+00, .1542E+01,
	773	& -.6903E-04,-.1081E-02,-.2849E-01,-.2557E+00,-.4515E+00,
	774	& -.6935E+00,-.9313E+00,-.1098E+01,-.1165E+01,-.1222E+01,
	775	& -.1239E+01,-.1285E+01,-.1353E+01,-.1432E+01,-.1606E+01,
	776	& -.1809E+01,-.1613E+01,-.1139E+01,-.1174E+01,-.1452E+01,
	777	& -.1075E+01,-.7739E+00,-.1122E+01,-.8658E+00, .8014E-02,
	778	& .2460E+00, .7393E+00, .1802E+01/
	779	DATA B1/-.4265E-01,-.8221E-01,-.9239E-01, .7328E+00, .1571E+01,
	780	& .2661E+01, .3805E+01, .4705E+01, .5192E+01, .5586E+01,
	781	& .5713E+01, .5844E+01, .5940E+01, .5939E+01, .5617E+01,
	782	& .6022E+01, .5165E+01, .2970E+01, .3196E+01, .4005E+01,
	783	& .1932E+01, .1164E+01, .2469E+01, .8862E+00,-.1572E+01,
	784	& -.1359E+01,-.2718E+01,-.5456E+01,
	785	& .2345E-03, .3756E-02, .1058E+00, .9941E+00, .1787E+01,
	786	& .2813E+01, .3909E+01, .4823E+01, .5400E+01, .5881E+01,
	787	& .6131E+01, .6373E+01, .6598E+01, .6769E+01, .6816E+01,
	788	& .7097E+01, .6380E+01, .4542E+01, .4294E+01, .5046E+01,
	789	& .3642E+01, .2361E+01, .3303E+01, .2337E+01,-.6819E+00,
	790	& -.1528E+01,-.3110E+01,-.6542E+01/
	791	DATA C1/ .5537E-01, .1094E+00, .1933E+00,-.3295E+00,-.9664E+00,
	792	& -.1836E+01,-.2791E+01,-.3579E+01,-.4033E+01,-.4384E+01,
	793	& -.4467E+01,-.4506E+01,-.4480E+01,-.4338E+01,-.3594E+01,
	794	& -.3592E+01,-.2857E+01,-.9067E+00,-.8569E+00,-.1483E+01,
	795	& .2426E+00, .1010E+01,-.2250E-01, .1272E+01, .3286E+01,
	796	& .3064E+01, .4112E+01, .6241E+01,
	797	& -.1499E-03,-.2538E-02,-.7950E-01,-.7859E+00,-.1437E+01,
	798	& -.2302E+01,-.3255E+01,-.4084E+01,-.4636E+01,-.5081E+01,
	799	& -.5287E+01,-.5434E+01,-.5522E+01,-.5524E+01,-.5105E+01,
	800	& -.5006E+01,-.4323E+01,-.2649E+01,-.2219E+01,-.2726E+01,
	801	& -.1514E+01,-.3276E+00,-.9927E+00,-.1745E+00, .2306E+01,
	802	& .2984E+01, .4190E+01, .6871E+01/
	803	DATA A2/ .8210E-03, .6438E-03,-.2567E-01,-.2701E+00,-.4978E+00,
	804	& -.8074E+00,-.1155E+01,-.1458E+01,-.1649E+01,-.1731E+01,
	805	& -.1683E+01,-.1675E+01,-.1646E+01,-.1544E+01,-.1234E+01,
	806	& -.1855E+01,-.1307E+01,-.2446E+00,-.5417E+00,-.1624E+01,
	807	& -.2091E+00, .2351E+00,-.1552E+01,-.5738E+00, .1313E+01,
	808	& .4078E-01, .1116E+01, .5153E+01,
	809	& -.7454E-04,-.1152E-02,-.2981E-01,-.2758E+00,-.5027E+00,
	810	& -.8111E+00,-.1158E+01,-.1463E+01,-.1658E+01,-.1747E+01,
	811	& -.1709E+01,-.1708E+01,-.1687E+01,-.1596E+01,-.1308E+01,
	812	& -.1910E+01,-.1371E+01,-.3399E+00,-.6009E+00,-.1651E+01,
	813	& -.3374E+00, .1413E+00,-.1566E+01,-.7407E+00, .1151E+01,
	814	& .9063E-01, .1134E+01, .5181E+01/
	815	DATA B2/-.4020E-02,-.4571E-02, .9034E-01, .1031E+01, .1929E+01,
	816	& .3185E+01, .4662E+01, .6060E+01, .7090E+01, .7687E+01,
	817	& .7703E+01, .7758E+01, .7644E+01, .7176E+01, .5232E+01,
	818	& .6690E+01, .4972E+01, .1198E+01, .1616E+01, .5180E+01,
	819	& .4740E+00,-.1502E+01, .4302E+01, .9918E+00,-.5496E+01,
	820	& -.1222E+01,-.4838E+01,-.1792E+02,
	821	& .2569E-03, .4046E-02, .1111E+00, .1065E+01, .1964E+01,
	822	& .3219E+01, .4696E+01, .6101E+01, .7148E+01, .7768E+01,
	823	& .7824E+01, .7907E+01, .7820E+01, .7389E+01, .5527E+01,
	824	& .6926E+01, .5219E+01, .1545E+01, .1852E+01, .5278E+01,
	825	& .8977E+00,-.1184E+01, .4303E+01, .1494E+01,-.5022E+01,
	826	& -.1536E+01,-.4991E+01,-.1814E+02/
	827	DATA C2/ .5370E-02, .8378E-02,-.5663E-01,-.7911E+00,-.1516E+01,
	828	& -.2549E+01,-.3786E+01,-.4980E+01,-.5871E+01,-.6372E+01,
	829	& -.6307E+01,-.6228E+01,-.5961E+01,-.5361E+01,-.3045E+01,
	830	& -.3777E+01,-.2316E+01, .8933E+00, .8813E+00,-.1974E+01,
	831	& .1810E+01, .3680E+01,-.9658E+00, .1690E+01, .7053E+01,
	832	& .3513E+01, .6392E+01, .1679E+02,
	833	& -.1729E-03,-.2839E-02,-.8489E-01,-.8456E+00,-.1578E+01,
	834	& -.2617E+01,-.3861E+01,-.5067E+01,-.5977E+01,-.6502E+01,
	835	& -.6476E+01,-.6425E+01,-.6184E+01,-.5615E+01,-.3369E+01,
	836	& -.4053E+01,-.2592E+01, .5416E+00, .6159E+00,-.2112E+01,
	837	& .1417E+01, .3365E+01,-.9937E+00, .1261E+01, .6647E+01,
	838	& .3797E+01, .6509E+01, .1696E+02/
	839	DATA itime /0/
	840	data xldo/0.1,0.5,1.0,2.0,4.0,8.0/
	841	data xx1/0.66,0.70,0.60,0.50,0.50,0.59/
	842	data xx2/0.66,0.70,0.81,0.60,0.59,0.62/
	843	data yy1/0.66,0.70,0.52,0.44,0.59,0.59/
	844	data yy2/0.66,0.70,0.74,0.55,0.62,0.62/
	845	data zz1/0.66,0.70,0.42,0.48,0.55,0.59/
	846	data zz2/0.66,0.70,0.65,0.55,0.63,0.62/
	847
	848	save ifirst
	849
	850	if (ifirst.eq.0) then
	851	print*,' IFIRST', ifirst
	852	do i=1,nech
	853	xech(i)=xech(i)*1.7
	854	enddo
	855	ifirst=1
	856	endif
	857
	858	** 1: compute the exponent of the law N^a LOG N
	859	** with the index n and the size parameter x
	860	****************************************************
	861	xexp5=0.66666
	862	xexp6=0.66666
	863	****** INTERPOLATION WITH THE VALUE OF REAL REFR. INDEX
	864	do i=1,nex
	865	if(xn.lt.1.7) then
	866	ww1(i)=(yy1(i)-xx1(i))/.3*(xn-1.4)+xx1(i)
	867	ww2(i)=(yy2(i)-xx2(i))/.3*(xn-1.4)+xx2(i)
	868	endif
	869	if(xn.ge.1.7) then
	870	ww1(i)=(zz1(i)-yy1(i))/.3*(xn-1.7)+yy1(i)
	871	ww2(i)=(zz2(i)-yy2(i))/.3*(xn-1.7)+yy2(i)
	872	endif
	873	enddo
	874	****** INTERPOLATION WITH THE SIZE PARAMETER
	875	******** XEXP5 AND XEXP6 ARE THE EXPONENT TO USE IN N^a LOGN LAW
	876	do i=2,nex
	877	if(x.lt.xldo(i).and.x.ge.xldo(i-1)) then
	878	rap=xldo(i)-xldo(i-1)
	879	xexp5=(ww1(i)-ww1(i-1))/rap*(x-xldo(i-1))+ww1(i-1)
	880	xexp6=(ww2(i)-ww2(i-1))/rap*(x-xldo(i-1))+ww2(i-1)
	881	endif
	882	enddo
	883	if(x.ge.xldo(nex)) then
	884	xexp5=ww1(nex)
	885	xexp6=ww2(nex)
	886	endif
	887
	888	xxn=xn*x
	889
	890	** location in the defined range of n
	891	********************************************
	892	** Out of the range
	893	xni=xn
	894	dxn=0.
	895	if (xn.lt.1.4) then
	896	dxn=(xn-1.4)
	897	xn=1.4
	898	endif
	899	if (xn.gt.2.0) then
	900	dxn=(xn-2.0)
	901	xn=2.0
	902	endif
	903
	904	** location in the defined ranges of n*X
	905	********************************************
	906	** Out of the range
	907	if (xxn.ge.xech(nech)) xxn=xech(nech)*.99
	908	** Rayleigh range
	909	if (xxn.lt.xech(1)) xxn=xech(1)
	910	do i=1,nech
	911	if(xech(i).le.xxn) j=i
	912	enddo
	913	** Location inside the slab j to j+1
	914	xechj=xech(j)
	915	xechjp1=xech(j+1)
	916	jindex=j
	917	xrap=(xxn-xech(j))/(xech(j+1)-xech(j))
	918	xrapl=(alog10(xxn)-alog10(xech(j)))
	919	& /(alog10(xech(j+1))-alog10(xech(j)))
	920	if(xxn.gt.1.7) xrapl=xrap
	921	***************************************************
	922	Calculation of the (parabolic) coefficients
	923	***************************************************
	924	xki=xk
	925	if (xk.lt.1.e-2) xki=1.e-2
	926	if (xk.gt.1.) xki=1.
	927	rki=-alog10(xki)
	928	**** Computation of the parabolic coefficients
	929	** For index j
	930	* f(x0+dx)=f(x0)+[df/dx](x0)*dx
	931	* with f=ax^2+bx+c ---> f'=[f-c+ax^2]/x=2ax+b
	932	*
	933	if(rki.gt.2.) then
	934	cjd =A2(j) xn2.+B2(j) xn+C2(j)
	935	cjd_=A2(j+nech)xn2.+B2(j+nech)xn+C2(j+nech)
	936	cju =A2(j) xn2.+B2(j) xn+C2(j)
	937	cju_=A2(j+nech)xn2.+B2(j+nech)xn+C2(j+nech)
	938	* neutre si dxn=0 :
	939	cjd =cjd +(cjd -C2(j) +A2(j) xn2.)/xndxn
	940	cjd_=cjd_+(cjd_-C2(j+nech)+A2(j+nech)xn2.)/xndxn
	941	cju =cju +(cju -C2(j) +A2(j) xn2.)/xndxn
	942	cju_=cju_+(cju_-C2(j+nech)+A2(j+nech)xn2.)/xndxn
	943	xku=1.e-2
	944	xkd=1.e-2
	945	endif
	946	if(rki.le.2.) then
	947	cjd =A1(j) xn2.+B1(j) xn+C1(j)
	948	cjd_=A1(j+nech)xn2.+B1(j+nech)xn+C1(j+nech)
	949	cju =A2(j) xn2.+B2(j) xn+C2(j)
	950	cju_=A2(j+nech)xn2.+B2(j+nech)xn+C2(j+nech)
	951	* neutre si dxn=0 :
	952	cjd =cjd +(cjd -C1(j) +A1(j) xn2.)/xndxn
	953	cjd_=cjd_+(cjd_-C1(j+nech)+A1(j+nech)xn2.)/xndxn
	954	cju =cju +(cju -C2(j) +A2(j) xn2.)/xndxn
	955	cju_=cju_+(cju_-C2(j+nech)+A2(j+nech)xn2.)/xndxn
	956	xku=1.e-2
	957	xkd=1.e-1
	958	endif
	959	if(rki.le.1.) then
	960	cju =A1(j) xn2.+B1(j) xn+C1(j)
	961	cju_=A1(j+nech)xn2.+B1(j+nech)xn+C1(j+nech)
	962	cjd =A0(j) xn2.+B0(j) xn+C0(j)
	963	cjd_=A0(j+nech)xn2.+B0(j+nech)xn+C0(j+nech)
	964	* neutre si dxn=0 :
	965	cjd =cjd +(cjd -C0(j) +A0(j) xn2.)/xndxn
	966	cjd_=cjd_+(cjd_-C0(j+nech)+A0(j+nech)xn2.)/xndxn
	967	cju =cju +(cju -C1(j) +A1(j) xn2.)/xndxn
	968	cju_=cju_+(cju_-C1(j+nech)+A1(j+nech)xn2.)/xndxn
	969	xku=1.e-1
	970	xkd=1.e-0
	971	endif
	972
	973
	974	** For index j+1
	975	if(rki.gt.2.) then
	976	cjdp1 =A2(j+1) xn2.+B2(j+1) xn+C2(j+1)
	977	cjdp1_=A2(j+1+nech)xn2.+B2(j+1+nech)xn+C2(j+1+nech)
	978	cjup1 =A2(j+1) xn2.+B2(j+1) xn+C2(j+1)
	979	cjup1_=A2(j+1+nech)xn2.+B2(j+1+nech)xn+C2(j+1+nech)
	980	* neutre si dxn=0 :
	981	cjdp1 =cjdp1 +(cjdp1 -C2(j+1) +A2(j+1) xn2.)/xndxn
	982	cjdp1_=cjdp1_+(cjdp1_-C2(j+1+nech)+A2(j+1+nech)xn2.)/xndxn
	983	cjup1 =cjup1 +(cjup1 -C2(j+1) +A2(j+1) xn2.)/xndxn
	984	cjup1_=cjup1_+(cjup1_-C2(j+1+nech)+A2(j+1+nech)xn2.)/xndxn
	985	xku=1.e-2
	986	xkd=1.e-2
	987	endif
	988	if(rki.le.2.) then
	989	cjdp1 =A1(j+1) xn2.+B1(j+1) xn+C1(j+1)
	990	cjdp1_=A1(j+1+nech)xn2.+B1(j+1+nech)xn+C1(j+1+nech)
	991	cjup1 =A2(j+1) xn2.+B2(j+1) xn+C2(j+1)
	992	cjup1_=A2(j+1+nech)xn2.+B2(j+1+nech)xn+C2(j+1+nech)
	993	* neutre si dxn=0 :
	994	cjdp1 =cjdp1 +(cjdp1 -C1(j+1) +A1(j+1) xn2.)/xndxn
	995	cjdp1_=cjdp1_+(cjdp1_-C1(j+1+nech)+A1(j+1+nech)xn2.)/xndxn
	996	cjup1 =cjup1 +(cjup1 -C2(j+1) +A2(j+1) xn2.)/xndxn
	997	cjup1_=cjup1_+(cjup1_-C2(j+1+nech)+A2(j+1+nech)xn2.)/xndxn
	998	xku=1.e-2
	999	xkd=1.e-1
	1000	endif
	1001	if(rki.le.1.) then
	1002	cjup1 =A1(j+1) xn2.+B1(j+1) xn+C1(j+1)
	1003	cjup1_=A1(j+1+nech)xn2.+B1(j+1+nech)xn+C1(j+1+nech)
	1004	cjdp1 =A0(j+1) xn2.+B0(j+1) xn+C0(j+1)
	1005	cjdp1_=A0(j+1+nech)xn2.+B0(j+1+nech)xn+C0(j+1+nech)
	1006	* neutre si dxn=0 :
	1007	cjdp1 =cjdp1 +(cjdp1 -C0(j+1) +A0(j+1) xn2.)/xndxn
	1008	cjdp1_=cjdp1_+(cjdp1_-C0(j+1+nech)+A0(j+1+nech)xn2.)/xndxn
	1009	cjup1 =cjup1 +(cjup1 -C1(j+1) +A1(j+1) xn2.)/xndxn
	1010	cjup1_=cjup1_+(cjup1_-C1(j+1+nech)+A1(j+1+nech)xn2.)/xndxn
	1011	xku=1.e-1
	1012	xkd=1.e-0
	1013	endif
	1014
	1015
	1016	** Computation of the monomer-factor
	1017
	1018	cjup1 =cjup1*20.
	1019	cjup1_=cjup1_*20.
	1020	cju =cju*20.
	1021	cju_ =cju_*20.
	1022	cjdp1 =cjdp1*20.
	1023	cjdp1_=cjdp1_*20.
	1024	cjd =cjd*20.
	1025	cjd_ =cjd_*20.
	1026	xn=xni ! restitution de xn
	1027
	1028	return
	1029	end
	1030
	1031	subroutine structure(NB,X,Z,STRUCT,DF)
	1032	implicit real (a-h,o-z)
	1033	c integer NB
	1034	real NB
	1035	real X,Z,D
	1036
	1037	D=DF
	1038	if (DF.eq.2) D=2.0001
	1039	if (z.eq.0.) z=1.e-4
	1040
	1041	STRUCT=1.
	1042	NOSTRUCT=0 ! if asymetry parameters are not needed
	1043	! just skip the computation: save your time!
	1044	if (NOSTRUCT.eq.1) goto 102
	1045	u0=5.
	1046	pi=3.1415926
	1047	* If convergence test in on (end of the loop):
	1048	xacc=1.e-3
	1049	* Else, computation is done once: accuracy is generally about 1%
	1050
	1051	* The structure factor is computed in order to evaluate the asymetry
	1052	* parameter (not for cross sections). We need to compute the integral
	1053	* of the following function:
	1054	*
	1055	* sin(2XZu)exp(-1/2u**2) for u between 0 and 5.
	1056	*
	1057	* where X,Z are provided through the subroutine calling
	1058	* A=4*pi (normalization factor for D=2 --> Botet et al., 1995)
	1059	*
	1060	* And STRUCT is given as:
	1061	*
	1062	* STRUCT=N[1+(N-1)2pi/(AXZ) INTEGRAL[sin(2XZu)exp(-1/2u*2)du]
	1063	*
	1064	* The number of oscillations for sin(2XZu) between 0 and U is:
	1065	* n=UXZ/pi.... let integer (simpson integration).
	1066
	1067	A=-5.026*D+22.618
	1068	A=A/(2.pi) !<---- here is A/(2PI)
	1069	nplt=6int(5.Z*X/3.1415926+1.)
	1070	* nplt=int(5.ZX/3.1415926+1.)
	1071	* This is the number of periods for the sinus in the range
	1072	* u E [0,5]. Integration is done with 6 points per period.
	1073	* Accuracy is about 1% on the final result STRUCT.
	1074	*
	1075	* The minimum value for nplt is set to 17 to do computation
	1076	* in the Rayleigh range ( when Z*X reaches 0)
	1077	STRUCT=1.e-10
	1078	101 if ((nplt/2)1..eq.nplt.5) nplt=nplt+1 !---> odd
	1079	if (nplt.lt.17) nplt=17
	1080	dint=u0/(nplt-1)
	1081	STRUCT_OLD=STRUCT
	1082	STRUCT=0.
	1083
	1084	*** EXTENDED SIMPSON INTEGRATION
	1085	iint=0
	1086	ucr=iint*dint
	1087	um1=sin(2.xzucr)exp(-.5ucrD)ucr**(D-2.)
	1088	STRUCT=STRUCT+um1
	1089	iint=nplt-1
	1090	ucr=iint*dint
	1091	um1=sin(2.xzucr)exp(-.5ucrD)ucr**(D-2.)
	1092	STRUCT=STRUCT+um1
	1093
	1094	do iint=1,nplt-2,2
	1095	ucr=iint*dint
	1096	um1=4.sin(2.xzucr)exp(-.5ucr*D)ucr**(D-2.)
	1097	STRUCT=STRUCT+um1
	1098	enddo
	1099
	1100	do iint=2,nplt-3,2
	1101	ucr=iint*dint
	1102	um1=2.sin(2.xzucr)exp(-.5ucr*D)ucr**(D-2.)
	1103	STRUCT=STRUCT+um1
	1104	enddo
	1105	STRUCT=dint/3.*STRUCT
	1106	ERR=abs(STRUCT_OLD-STRUCT)/STRUCT
	1107	nplt=int(nplt*2)
	1108	C UNCOMMENT THE IF STATEMENT TO
	1109	c SET ON THE CONVERGENCE TEST:
	1110	c if (ERR.gt.xacc) GOTO 101
	1111	STRUCT=(STRUCT/(xza)(NB-1)+1.)NB
	1112
	1113	102 continue
	1114	return
	1115	end
	1116
	1117

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format