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

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

fix very old (!) bug in psvfcj.for
implement neutron scattering lengths in GSASII including the dozen anomalous scattering isotopes
isotope choice is in General
so GSASII will now do neutron CW Rietveld refinements
some cleanup of the constraints GUI
remove varyList from GSASIImapvars.py globals

File size: 10.1 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)
129        tmp2 = 1.0 + ApB2
130        tmp = SQRT(tmp2)*cos2THETA
131        if (abs(tmp) .le. 1.0) then
132          Emin = Acosd(tmp)
133          tmp1 = tmp2*(1.0 - tmp2*(1.0-sin2THETA2))
134        else
135            print *,'tmp > 1.0'
136          tmp1 = 0.0
137          if (tmp .gt. 0.0) then
138            Emin = 0.0
139          else
140            Emin = 180.0
141          endif
142        endif
143        if (tmp1 .gt. 0 .and. abs(tmp) .le. 1.0) then
144          dEmindA = -ApB*cos2THETA/SQRT(tmp1)
145        ELSE
146          dEmindA = 0.0
147        ENDIF
148
149c       
150c Compute Gaussian Quadrature interval
151C Note that the number of points must be large enough so that the
152c interval between 2phi(min) and 2theta must be divided into steps no
153c larger than 0.005 degrees.  LWF  8/10/95
154c       
155c Determine which Gauss-Legendre Table to use
156c       
157        arraynum = 1
158c       
159c Calculate desired number of intervals
160c
161        tmp = abs(TTHETAD - emin)               ! tmp in degrees
162        IF ( GAM.LE.0.0 ) THEN
163          K = INT(TMP*200.0)                    !RVD formulation
164        ELSE
165          k = int(300.0*tmp/GAM)                !New formulation of May 2004
166        END IF
167c
168C Find the next largest set
169c
170        do while (arraynum .lt. numtab .and. k.gt.nterms(arraynum))
171          arraynum = arraynum + 1
172        enddo
173        ngt = nterms(arraynum)
174c
175C calculate the terms, if they have not been used before
176c
177        if (.not. calclfg(arraynum)) then
178          calclfg(arraynum) = .true.
179          call gauleg(-1.,1.,xpt,wpt,ngt)
180          it = fstterm(arraynum)-ngt/2
181c
182C copy the ngt/2 terms from our working array to the stored array
183c
184          do k=ngt/2+1,ngt
185            xp(k+it) = xpt(k)
186            wp(k+it) = wpt(k)
187          enddo
188        endif
189        sumWG = 0.
190        sumWRG = 0.
191        sumWdGdA = 0.
192        sumWRdGdA = 0.
193        sumWGdRd2t = 0.
194        sumWGdRdsig = 0.
195        sumWGdRdgam = 0.
196        sumWGdRdA = 0.
197        sumWGdRdB = 0.
198c       
199c Compute Convolution integral for 2phi(min) <= delta <= 2theta
200c       
201        it = fstterm(arraynum)-ngt/2
202        do k=ngt/2+1,ngt
203          delta = emin + (TTHETAD - emin) * xp(k+it) ! Delta in degrees
204          CALL PSVOIGT(DTT+TTHETA-delta*100.,SIG,GAM,R,dRdT,dRdS,dRdG)
205
206          dDELTAdA = (1. - xp(k+it))*dEmindA ! N. B. symm w/r A,B
207          sinDELTA = sind(Delta)
208          cosDELTA = cosd(Delta)
209          RcosDELTA = 1. / cosDELTA
210          tanDELTA = tand(Delta)
211          cosDELTA2 = cosDELTA*cosDELTA     
212          if ( ttheta.le.4500.0 .or. ttheta.gt.13500.0 ) then
213            tmp = sin2THETA2-sinDelta**2
214          else
215            tmp = cosDELTA2-cos2THETA2
216          end if
217          if (tmp .gt. 0) then
218            tmp1 = 1.0/SQRT(tmp)
219            F = abs(cos2THETA) * tmp1
220            dFdA = cosDELTA*cos2THETA*sinDELTA*dDELTAdA*tmp1**3
221          else
222            F = 1.0
223            dFdA = 0.0
224          endif
225c       
226c Calculate G(Delta,2theta) [G = W /(h cos(delta) ] [ FCJ eq. 7(a) and 7(b) ]
227c
228          if ( ttheta.gt.9000.0 ) then                  !are you kidding!! this worked
229!          if(abs(delta-emin) .gt. abs(einfl-emin))then
230              G = 2.0*A*F*RcosDELTA
231              dGdA = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
232          else
233            G = (-1.0 + ApB*F) * RcosDELTA
234            dGdA = RcosDELTA*(F - tanDELTA*dDELTAdA
235     1             + ApB*F*tanDELTA*dDELTAdA + ApB*dFdA)
236          endif
237
238          WG = wp(k+it) * G
239          sumWG = sumWG + WG
240          sumWRG = sumWRG + WG * R
241          sumWdGdA = sumWdGdA + wp(k+it) * dGdA
242          sumWRdGdA = sumWRdGdA + wp(k+it) * R * dGdA
243          sumWGdRd2t = sumWGdRd2t + WG * dRdT ! N.B. 1/centidegrees
244          sumWGdRdsig = sumWGdRdsig + WG * dRdS
245          sumWGdRdgam = sumWGdRdgam + WG * dRdG
246          sumWGdRdA = sumWGdRdA + WG * dRdT * dDELTAdA ! N. B. symm w/r A,B
247        enddo
248        RsumWG = 1.0/sumWG
249        RsumWG2 = RsumWG * RsumWG
250        PRFUNC = sumWRG * RsumWG
251
252        dydA = (-(sumWRG*sumWdGdA) +
253     1    sumWG*(sumWRdGdA- 100.0 * todeg * sumWGdRdA)) * RsumWG2
254        sigpart = sumWGdRdsig * RsumWG
255        gampart = sumWGdRdgam * RsumWG
256        DPRDT = -SumWGdRd2T * RsumWG
257      else
258c
259C no asymmetry -- nice and simple!
260c
261        CALL PSVOIGT(DTT,SIG,GAM,R,dRdT,dRdS,dRdG)
262        PRFUNC = R
263        dydA = 0.004 * sign(1.0,TTHETA - DTT)
264        sigpart = dRdS
265        gampart = dRdG
266        DPRDT = -dRdT
267      END IF
268      SHLPART = DYDA
269     
270      RETURN
271      END
272     
Note: See TracBrowser for help on using the repository browser.