source: trunk/fsource/powsubs/psvfcj.for @ 209

Last change on this file since 209 was 209, checked in by vondreele, 13 years ago

fortran fixes - new self contained libraries
had to change all INCLUDE copyright lines.
File size: 11.4 KB
Line 
1      SUBROUTINE PSVFCJ(DTT,TTHETA,SL,HL,SIG,GAM,PRFUNC,DPRDT,SLPART,
2     1  HLPART,SIGPART,GAMPART)
3
4!PURPOSE: Compute function & derivatives for Pseudovoigt profile
5!   [W.I.F.David (1986), J. Appl. Cryst. 19, 63-64 &
6!    P. Thompson, D.E. Cox & J.B. Hastings (1987) J. Appl. Cryst.,20,79-83.]
7! Finger-Cox-Jephcoat (FCJ94) asymmetry correction
8!   [L.W. Finger, D.E. Cox & A.P. Jephcoat (1994) J. Appl. Cryst.,27,892-900.]
9! coded 11/95 by B. H. Toby (NIST). revised version
10! parameterized as asym1=S/L asym2=H/L
11
12
13      INCLUDE       'INCLDS/COPYRIGT.FOR' 
14
15!CALLING ARGUMENTS:
16
17      REAL*4        DTT                 !delta 2-theta in centidegrees                 
18      REAL*4        TTHETA              !2-theta in centidegrees             
19      REAL*4        SL,HL               !S/L & H/L               
20      REAL*4        SIG,GAM             
21      REAL*4        PRFUNC             
22      REAL*4        DPRDT               
23      REAL*4        SLPART,HLPART       
24      REAL*4        SIGPART,GAMPART     
25
26!INCLUDE STATEMENTS:
27      real*4 sind,cosd,tand,acosd
28
29!LOCAL VARIABLES:
30
31      REAL*4        R                   ! pseudo-voight intensity
32      REAL*4        DRDT                ! deriv R w/r theta
33      REAL*4        DRDS                ! deriv R w/r sig
34      REAL*4        DRDG                ! deriv R w/r gam
35      REAL*4        F                   
36      REAL*4        G                   
37      REAL*4        DFDA               
38      REAL*4        DGDA               
39      REAL*4        DGDB               
40      REAL*4        DYDA               
41      REAL*4        DYDB               
42      REAL*4        SIN2THETA2          ! sin(2theta)**2
43      REAL*4        COS2THETA           ! cos(2theta)
44      REAL*4        SIN2THETA           ! sin(2THETA)
45      REAL*4        SINDELTA            ! sin(Delta)
46      REAL*4        COSDELTA            ! cos(Delta)
47      REAL*4        RCOSDELTA           ! 1/cos(Delta)
48      REAL*4        TANDELTA            ! tan(Delta)
49      REAL*4        COSDELTA2           ! cos(Delta)**2
50      REAL*4        A                   ! asym1 [coff(7)]
51      REAL*4        B                   ! asym2 [coff(8)]
52      REAL*4        APB                 ! (A+B)
53      REAL*4        AMB                 ! (A-B)
54      REAL*4        APB2                ! (A+B)**2
55      REAL*4        TTHETAD             ! Two Theta in degrees
56
57! Intermediate variables
58
59      REAL*4        RSUMWG2             !      1.0/(sum of w G)**2
60      REAL*4        SUMWG               !      sum of w G
61      REAL*4        WG                  !      w G
62      REAL*4        RSUMWG              !      1.0/sum of w G
63      REAL*4        SUMWRG              !      sum of w G
64      REAL*4        SUMWDGDA            !      sum of w dGdA
65      REAL*4        SUMWRDGDA           !      sum of w R dGdA
66      REAL*4        SUMWDGDB            !      sum of w dGdB
67      REAL*4        SUMWRDGDB           !      sum of w R dGdB
68      REAL*4        SUMWGDRD2T          !      sum of w G dRd(2theta)
69      REAL*4        SUMWGDRDSIG         !      sum of w G dRdp(n)
70      REAL*4        SUMWGDRDGAM         !      sum of w G dRdp(n)
71      REAL*4        SUMWGDRDA           
72      REAL*4        SUMWGDRDB           
73      REAL*4        EMIN                ! 2phi minimum
74      REAL*4        EINFL               ! 2phi of inflection point
75      REAL*4        DEMINDA             ! Derivative of Emin wrt A
76      REAL*4        DELTA               ! Angle of integration for convolution
77      REAL*4        DDELTADA           
78      REAL*4        TMP,TMP1,TMP2       ! intermediates
79      INTEGER*4     I,K,IT              ! Miscellaneous loop variables
80c       
81c       Local Variables for Gaussian Integration
82c       
83      INTEGER*4     NGT                 !number of terms in Gaussian quadrature
84      INTEGER*4     NUMTAB              ! number of pre-computed Gaussian tables
85      parameter      (numtab=34)
86      INTEGER*4     NTERMS(NUMTAB)      ! number of terms in each table - must be even
87      INTEGER*4     FSTTERM(NUMTAB)     ! location of 1st term: N.B. N/2 terms
88      LOGICAL*4     CALCLFG(NUMTAB)     ! true if table has previously been calculated
89      INTEGER*4     ARRAYNUM            ! number of selected array
90      INTEGER*4     ARRAYSIZE           ! size of complete array
91      parameter      (arraysize=1670)
92      REAL*4        XP(ARRAYSIZE)       !Gaussian abscissas
93      REAL*4        WP(ARRAYSIZE)       !Gaussian weights
94      REAL*4        XPT(400)            !temporary Gaussian abscissas
95      REAL*4        WPT(400)            !temporary Gaussian weights
96      REAL*4        STOFW               
97      PARAMETER (STOFW=2.35482005)
98      REAL*4        TODEG               
99      PARAMETER (todeg=57.2957795)
100      save      calclfg,xp,wp          !Values to be saved across calls
101
102!FUNCTION DEFINITIONS:
103
104!DATA STATEMENTS:
105
106      DATA      NTERMS                ! number of terms in each table - must be even
107     1  /2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,40,50,60,70,80,90     ! Numbers of terms changed May 2004
108     2   ,100,110,120,140,160,180,200,220,240,260,280,300,400/
109C       note that nterms determines both arraysize and fstterm
110      DATA      FSTTERM                ! loc. of 1st term: N.B. N/2 terms are saved
111     1  /0,1,3,6,10,15,21,28,36,45,55,66,78,91,105,120,140,165,195,230, !FSTERM(1) should be 0 - indexing starts at 1+this number!
112     2  270,315,365,420,480,550,630,720,820,930,1050,1180,1320,1470/
113      DATA      calclfg/numtab*.false./ ! true if table entry has been calculated
114
115!CODE:
116
117c
118C f(2theta) intermediates
119c
120      TTHETAD = 0.01 * TTHETA
121      SIN2THETA = SIND(TTHETAD)
122      cos2THETA = COSD(TTHETAD)
123      sin2theta2 = sin2THETA * sin2THETA
124      cos2THETA2 = cos2THETA * cos2THETA
125c
126C Asymmetry terms
127c
128      A = SL            ! A = S/L in FCJ
129      B = HL            ! B = H/L in FCJ
130      ApB = A+B
131      AmB = A-B
132      ApB2 = ApB*ApB
133c
134C handle the case where there is asymmetry
135c
136      IF (A .ne. 0.0 .or. B .ne. 0.0) then
137        Einfl = Acosd(SQRT(1.0 + AmB**2)*cos2THETA) ! 2phi(infl) FCJ eq 5 (degrees)
138        tmp2 = 1.0 + ApB2
139        tmp = SQRT(tmp2)*cos2THETA
140c
141C Treat case where A or B is zero - Set Einfl = 2theta
142c
143        if (A.eq.0.0 .or. B .eq. 0.0)Einfl = Acosd(cos2THETA)
144        if (abs(tmp) .le. 1.0) then
145          Emin = Acosd(tmp)      ! 2phi(min) FCJ eq 4 (degrees)
146          tmp1 = tmp2*(1.0 - tmp2*(1.0-sin2THETA2))
147        else
148          tmp1 = 0.0
149          if (tmp .gt. 0.0) then
150            Emin = 0.0
151          else
152            Emin = 180.0
153          endif
154        endif
155        if (tmp1 .gt. 0 .and. abs(tmp) .le. 1.0) then
156          dEmindA = -ApB*cos2THETA/SQRT(tmp1) ! N. B. symm w/r A,B
157        ELSE
158          dEmindA = 0.0
159        ENDIF
160
161c       
162c Compute Gaussian Quadrature interval
163C Note that the number of points must be large enough so that the
164c interval between 2phi(min) and 2theta must be divided into steps no
165c larger than 0.005 degrees.  LWF  8/10/95
166c       
167c Determine which Gauss-Legendre Table to use
168c       
169        arraynum = 1
170c       
171c Calculate desired number of intervals
172c
173        tmp = abs(TTHETAD - emin)               ! tmp in degrees
174        IF ( GAM.LE.0.0 ) THEN
175          K = INT(TMP*200.0)                    !RVD formulation
176        ELSE
177          k = int(300.0*tmp/GAM)                !New formulation of May 2004
178        END IF
179c
180C Find the next largest set
181c
182        do while (arraynum .lt. numtab .and. k.gt.nterms(arraynum))
183          arraynum = arraynum + 1
184        enddo
185        ngt = nterms(arraynum)
186c
187C calculate the terms, if they have not been used before
188c
189        if (.not. calclfg(arraynum)) then
190          calclfg(arraynum) = .true.
191          call gauleg(-1.,1.,xpt,wpt,ngt)
192          it = fstterm(arraynum)-ngt/2
193c
194C copy the ngt/2 terms from our working array to the stored array
195c
196          do k=ngt/2+1,ngt
197            xp(k+it) = xpt(k)
198            wp(k+it) = wpt(k)
199          enddo
200        endif
201        sumWG = 0.
202        sumWRG = 0.
203        sumWdGdA = 0.
204        sumWRdGdA = 0.
205        sumWdGdB = 0.
206        sumWRdGdB = 0.
207        sumWGdRd2t = 0.
208        sumWGdRdsig = 0.
209        sumWGdRdgam = 0.
210        sumWGdRdA = 0.
211        sumWGdRdB = 0.
212c       
213c Compute Convolution integral for 2phi(min) <= delta <= 2theta
214c       
215        it = fstterm(arraynum)-ngt/2
216        do k=ngt/2+1,ngt
217          delta = emin + (TTHETAD - emin) * xp(k+it) ! Delta in degrees
218          CALL PSVOIGT(DTT+TTHETA-delta*100.,SIG,GAM,R,dRdT,dRdS,dRdG)
219
220          dDELTAdA = (1. - xp(k+it))*dEmindA ! N. B. symm w/r A,B
221          sinDELTA = sind(Delta)
222          cosDELTA = cosd(Delta)
223          RcosDELTA = 1. / cosDELTA
224          tanDELTA = tand(Delta)
225          cosDELTA2 = cosDELTA*cosDELTA     
226          tmp1 = cosDELTA2 - cos2THETA2
227          tmp2 = sin2THETA2 - sinDelta * sinDELTA
228          tmp = tmp2
229          if ( ttheta.gt.4500.0 ) tmp = tmp1
230          if (tmp .gt. 0) then
231            tmp1 = 1.0/SQRT(tmp)
232            F = abs(cos2THETA) * tmp1
233            dFdA = cosDELTA*cos2THETA*sinDELTA*dDELTAdA * 
234     1           (tmp1*tmp1*tmp1)
235          else
236            F = 1.0
237            dFdA = 0.0
238          endif
239c       
240c Calculate G(Delta,2theta) [G = W /(h cos(delta) ] [ FCJ eq. 7(a) and 7(b) ]
241c       
242          if(abs(delta-emin) .gt. abs(einfl-emin))then
243            if ( A.ge.B) then
244c
245C N.B. this is the only place where d()/dA <> d()/dB
246c
247              G = 2.0*B*F*RcosDELTA
248              dGdA = 2.0*B*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
249              dGdB = dGdA + 2.0*F*RcosDELTA
250            else
251              G = 2.0*A*F*RcosDELTA
252              dGdB = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
253              dGdA = dGdB + 2.0*F*RcosDELTA
254            endif
255          else                                            ! delta .le. einfl .or. min(A,B) .eq. 0
256            G = (-1.0 + ApB*F) * RcosDELTA
257            dGdA = RcosDELTA*(F - tanDELTA*dDELTAdA
258     1             + ApB*F*tanDELTA*dDELTAdA + ApB*dFdA)
259            dGdB = dGdA
260          endif
261
262          WG = wp(k+it) * G
263          sumWG = sumWG + WG
264          sumWRG = sumWRG + WG * R
265          sumWdGdA = sumWdGdA + wp(k+it) * dGdA
266          sumWRdGdA = sumWRdGdA + wp(k+it) * R * dGdA
267          sumWdGdB = sumWdGdB + wp(k+it) * dGdB
268          sumWRdGdB = sumWRdGdB + wp(k+it) * R * dGdB
269          sumWGdRd2t = sumWGdRd2t + WG * dRdT ! N.B. 1/centidegrees
270          sumWGdRdsig = sumWGdRdsig + WG * dRdS
271          sumWGdRdgam = sumWGdRdgam + WG * dRdG
272          sumWGdRdA = sumWGdRdA + WG * dRdT * dDELTAdA ! N. B. symm w/r A,B
273        enddo
274        RsumWG = 1.0/sumWG
275        RsumWG2 = RsumWG * RsumWG
276        PRFUNC = sumWRG * RsumWG
277
278        dydA = (-(sumWRG*sumWdGdA) +
279     1    sumWG*(sumWRdGdA- 100.0 * todeg * sumWGdRdA)) * RsumWG2
280        dydB = (-(sumWRG*sumWdGdB) +
281     1    sumWG*(sumWRdGdB- 100.0 * todeg * sumWGdRdA)) * RsumWG2
282        sigpart = sumWGdRdsig * RsumWG
283        gampart = sumWGdRdgam * RsumWG
284        DPRDT = -SumWGdRd2T * RsumWG
285      else
286c
287C no asymmetry -- nice and simple!
288c
289        CALL PSVOIGT(DTT,SIG,GAM,R,dRdT,dRdS,dRdG)
290        PRFUNC = R
291        dydA = 0.002 * sign(1.0,TTHETA - DTT)
292        dydB = dydA
293        sigpart = dRdS
294        gampart = dRdG
295        DPRDT = -dRdT
296      END IF
297      SLPART = DYDA
298      HLPART = DYDB
299     
300      RETURN
301      END
302     
Note: See TracBrowser for help on using the repository browser.