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

Last change on this file since 353 was 353, checked in by vondreele, 10 years ago

modify psvfcj.for to use sh/l only
add derivative routine to pypowder.for
GSASIIgrid.py - new LS controls for derivative type
GSASIIlattice.py - work on docs as per P Jemian
GSASIIpwd.py - major mods to use FCJpsvoight fortran code & derivatives
GSASIIpwdGUI.py - mod to use LS controls

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