source: trunk/LMDZ.TITAN/libf/phytitan/cfffv11.F @ 537

Last change on this file since 537 was 175, checked in by slebonnois, 13 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
Line 
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
36c      refrel=(xn,xk)
37       refrel=cmplx(xn,xk)
38       nindex=int(xn*100.)               
39       x=2.*3.14159265*rad/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
88c      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
101c      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
115c      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.*pi*x**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
158c      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
240c         refrel=(xn,1.)
241          refrel=cmplx(xn,1.)
242          call intmie(x,refrel,nang,s1,s2,qextmd_,qscamd_)
243          qabsmd_=qextmd_-qscamd_
244
245c         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*(NC*1.)**.3333333
263          xsphe=x*(NC*1.)**.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(qs*NC*(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_=alphaS*x*(NB*1.)**(1./DFS)
351         call structure(NB,xx_,z,str,DFS)
352         if (z.eq.0.) str=(NB*1.)**2.
353         strS(j)=str
354  367   continue
355
356
357       do 365 j=1,2*nang-2,2
358
359        g1=2*dang/6*(s11(j)*sin(theta(j))*strS(j)
360     &              +s11(j+2)*sin(theta(j+2))*strS(j+2)
361     &            +4*s11(j+1)*sin(theta(j+1))*strS(j+1))
362     &                          *2*3.141592+g1
363
364        g2=2*dang/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     &      +4*s11(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=(NC*1.)**xexp5
394        rho_=(NB*1.)**xexp5
395
396* EQ 10 RESULT FOR ABSORPTION:
397*-----------------------------
398
399        sabs1=(sabs*xkf-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=(NTA*1.)**xexp5
410
411* COMPUTE sabs1 FOR NTA TO START THE POWER LAW:
412        sabs1=(sabs*xkf-qabs_)/alog10(20.)*alog10(NTA*1.)
413     &  *rhox/rho+qabs_
414
415* COMPUTE xexp1 EXPONENT FACTOR:
416                        xexp1=.85
417        if(x*xn.lt.2.2) xexp1=(.85-.96)*1.42*(x*xn-1.5)+.96
418        if(x*xn.lt.1.5) xexp1=.96
419
420* COMPUTE sabs1 VALUE FOR LARGE X:
421        sabs1=sabs1*(NB*1./(NTA*1.))**xexp1
422       endif
423
424* HERE SABS1 MEDIUM RANGE IS KNOWN
425
426        rho=(NC*1.)**xexp6
427        rho_=(NB*1.)**xexp6
428
429
430* EQ 10 RESULT FOR SCATTERING:
431*-----------------------------
432
433        ssca1=(ssca-qsca_)/alog10(20.)*alog10(NB*1.)
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=(NTS*1.)**xexp6
443
444* COMPUTE sabs1 FOR NTS TO START THE POWER LAW:
445        ssca1=(ssca-qsca_)/alog10(20.)*alog10(NTS*1.)
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*(NB*1./(NTS*1.))**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(NB*1.)+alog10(qsca_)
469        sext2=(alog10(sext)-alog10(qext_))/alog10(20.)
470     &            *alog10(NB*1.)+alog10(qext_)
471        sabs2=(alog10(sabs*xkf)-alog10(qabs_))/alog10(20.)
472     &            *alog10(NB*1.)+alog10(qabs_)
473
474
475
476************************************************
477* STEP # 7
478************************************************
479* THE CROSS OVER BETWEEN RAYLEYGH AND MEDIUM RANGE
480************************************************
481
482        rho_=(NB*1.)**.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=.5*alpha*x*(NB*1.)**.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.**(f*ssca1+g*ssca2)
517           fext=10.**(f*sext1+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=.5*alpha*x*(NB*1.)**.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.**(f*sabs1+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)*2*3.1415936*strS(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
583c new fashion...
584
585        fabs=(qext_-qsca_)*NB/(NB*1.)**(2./3.)*scoef
586        fext=fsca+fabs
587 
588
589        ENDIF
590
591        return
592       end
593
594c------------------------------------------------------------
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+4*x**.3333+2.
613       nstop=xstop
614       ymod=cabs(y)
615       nmx=amax1(xstop,ymod)+15
616c      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))
621c      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
637c        s_(j)=(0.,0.)
638c        s1(j)=(0.,0.)
639c777     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
649c      xi0=(apsi0,-chi0)
650c      xi1=(apsi1,-chi1)
651       xi0=cmplx(apsi0,-chi0)
652       xi1=cmplx(apsi1,-chi1)
653       qsca=0.
654       qsca_=0.
655       n=1
656c      *************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
664c      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=(refrel*d(n)+rn/x)*apsi-apsi1
669       bn=bn/((refrel*d(n)+rn/x)*xi-xi1)
670       qsca=qsca+(2*rn+1.)*(cabs(an)*cabs(an)+cabs(bn)*cabs(bn))
671 
672c     ***************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)=rn*amu(j)*pi(j)-(rn+1.)*pi0(j)  ! fonction de legendre
677          s1(j)=s1(j)+fn*(an*pi(j)+bn*tau(j))
678          s2(j)=s2(j)+fn*(an*tau(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*(an*pi(j)*p+bn*tau(j)*t)
683         s2(jj)=s2(jj)+fn*(an*tau(j)*t+bn*pi(j)*p)
684 789   continue
685       psi0=psi1
686       psi1=psi
687       apsi1=psi1           ! double prec=simple
688       chi0=chi1
689       chi1=chi
690c      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./(x*x))*qsca
700       qext=(4./(x*x))*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(2*nech), B0(2*nech), C0(2*nech), D0(2*nech)
713      real A1(2*nech), B1(2*nech), C1(2*nech), D1(2*nech)
714      real A2(2*nech), B2(2*nech), C2(2*nech), D2(2*nech)
715      real A3(2*nech), B3(2*nech), C3(2*nech), D3(2*nech)
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)     *xn**2.+B2(j)     *xn+C2(j)
935        cjd_=A2(j+nech)*xn**2.+B2(j+nech)*xn+C2(j+nech)
936        cju =A2(j)     *xn**2.+B2(j)     *xn+C2(j)
937        cju_=A2(j+nech)*xn**2.+B2(j+nech)*xn+C2(j+nech)
938* neutre si dxn=0  : 
939        cjd =cjd +(cjd -C2(j)     +A2(j)     *xn**2.)/xn*dxn
940        cjd_=cjd_+(cjd_-C2(j+nech)+A2(j+nech)*xn**2.)/xn*dxn
941        cju =cju +(cju -C2(j)     +A2(j)     *xn**2.)/xn*dxn
942        cju_=cju_+(cju_-C2(j+nech)+A2(j+nech)*xn**2.)/xn*dxn
943         xku=1.e-2
944         xkd=1.e-2
945       endif
946       if(rki.le.2.) then
947        cjd =A1(j)     *xn**2.+B1(j)     *xn+C1(j)
948        cjd_=A1(j+nech)*xn**2.+B1(j+nech)*xn+C1(j+nech)
949        cju =A2(j)     *xn**2.+B2(j)     *xn+C2(j)
950        cju_=A2(j+nech)*xn**2.+B2(j+nech)*xn+C2(j+nech)
951* neutre si dxn=0  :
952        cjd =cjd +(cjd -C1(j)     +A1(j)     *xn**2.)/xn*dxn
953        cjd_=cjd_+(cjd_-C1(j+nech)+A1(j+nech)*xn**2.)/xn*dxn
954        cju =cju +(cju -C2(j)     +A2(j)     *xn**2.)/xn*dxn
955        cju_=cju_+(cju_-C2(j+nech)+A2(j+nech)*xn**2.)/xn*dxn
956         xku=1.e-2
957         xkd=1.e-1
958       endif
959       if(rki.le.1.) then
960        cju =A1(j)     *xn**2.+B1(j)     *xn+C1(j)
961        cju_=A1(j+nech)*xn**2.+B1(j+nech)*xn+C1(j+nech)
962        cjd =A0(j)     *xn**2.+B0(j)     *xn+C0(j)
963        cjd_=A0(j+nech)*xn**2.+B0(j+nech)*xn+C0(j+nech)
964* neutre si dxn=0  :
965        cjd =cjd +(cjd -C0(j)     +A0(j)     *xn**2.)/xn*dxn
966        cjd_=cjd_+(cjd_-C0(j+nech)+A0(j+nech)*xn**2.)/xn*dxn
967        cju =cju +(cju -C1(j)     +A1(j)     *xn**2.)/xn*dxn
968        cju_=cju_+(cju_-C1(j+nech)+A1(j+nech)*xn**2.)/xn*dxn
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)     *xn**2.+B2(j+1)     *xn+C2(j+1)
977        cjdp1_=A2(j+1+nech)*xn**2.+B2(j+1+nech)*xn+C2(j+1+nech)
978        cjup1 =A2(j+1)     *xn**2.+B2(j+1)     *xn+C2(j+1)
979        cjup1_=A2(j+1+nech)*xn**2.+B2(j+1+nech)*xn+C2(j+1+nech)
980* neutre si dxn=0  :
981        cjdp1 =cjdp1 +(cjdp1 -C2(j+1)     +A2(j+1)     *xn**2.)/xn*dxn
982        cjdp1_=cjdp1_+(cjdp1_-C2(j+1+nech)+A2(j+1+nech)*xn**2.)/xn*dxn
983        cjup1 =cjup1 +(cjup1 -C2(j+1)     +A2(j+1)     *xn**2.)/xn*dxn
984        cjup1_=cjup1_+(cjup1_-C2(j+1+nech)+A2(j+1+nech)*xn**2.)/xn*dxn
985         xku=1.e-2
986         xkd=1.e-2
987       endif
988       if(rki.le.2.) then
989        cjdp1 =A1(j+1)     *xn**2.+B1(j+1)     *xn+C1(j+1)
990        cjdp1_=A1(j+1+nech)*xn**2.+B1(j+1+nech)*xn+C1(j+1+nech)
991        cjup1 =A2(j+1)     *xn**2.+B2(j+1)     *xn+C2(j+1)
992        cjup1_=A2(j+1+nech)*xn**2.+B2(j+1+nech)*xn+C2(j+1+nech)
993* neutre si dxn=0  :
994        cjdp1 =cjdp1 +(cjdp1 -C1(j+1)     +A1(j+1)     *xn**2.)/xn*dxn
995        cjdp1_=cjdp1_+(cjdp1_-C1(j+1+nech)+A1(j+1+nech)*xn**2.)/xn*dxn
996        cjup1 =cjup1 +(cjup1 -C2(j+1)     +A2(j+1)     *xn**2.)/xn*dxn
997        cjup1_=cjup1_+(cjup1_-C2(j+1+nech)+A2(j+1+nech)*xn**2.)/xn*dxn
998         xku=1.e-2
999         xkd=1.e-1
1000       endif
1001       if(rki.le.1.) then
1002        cjup1 =A1(j+1)     *xn**2.+B1(j+1)     *xn+C1(j+1)
1003        cjup1_=A1(j+1+nech)*xn**2.+B1(j+1+nech)*xn+C1(j+1+nech)
1004        cjdp1 =A0(j+1)     *xn**2.+B0(j+1)     *xn+C0(j+1)
1005        cjdp1_=A0(j+1+nech)*xn**2.+B0(j+1+nech)*xn+C0(j+1+nech)
1006* neutre si dxn=0  :
1007        cjdp1 =cjdp1 +(cjdp1 -C0(j+1)     +A0(j+1)     *xn**2.)/xn*dxn
1008        cjdp1_=cjdp1_+(cjdp1_-C0(j+1+nech)+A0(j+1+nech)*xn**2.)/xn*dxn
1009        cjup1 =cjup1 +(cjup1 -C1(j+1)     +A1(j+1)     *xn**2.)/xn*dxn
1010        cjup1_=cjup1_+(cjup1_-C1(j+1+nech)+A1(j+1+nech)*xn**2.)/xn*dxn
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)
1033c       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)*2*pi/(A*X*Z) 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/(2*PI)
1069         nplt=6*int(5.*Z*X/3.1415926+1.)
1070*        nplt=int(5.*Z*X/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.*x*z*ucr)*exp(-.5*ucr**D)*ucr**(D-2.)
1088          STRUCT=STRUCT+um1
1089         iint=nplt-1
1090         ucr=iint*dint
1091         um1=sin(2.*x*z*ucr)*exp(-.5*ucr**D)*ucr**(D-2.)
1092          STRUCT=STRUCT+um1
1093         
1094         do iint=1,nplt-2,2
1095         ucr=iint*dint
1096         um1=4.*sin(2.*x*z*ucr)*exp(-.5*ucr**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.*x*z*ucr)*exp(-.5*ucr**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)
1108C        UNCOMMENT THE IF STATEMENT TO
1109c        SET ON THE CONVERGENCE TEST:
1110c        if (ERR.gt.xacc) GOTO 101
1111         STRUCT=(STRUCT/(x*z*a)*(NB-1)+1.)*NB
1112
1113 102     continue       
1114         return
1115         end
1116
1117     
Note: See TracBrowser for help on using the repository browser.