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

Last change on this file since 2871 was 2871, checked in by vondreele, 6 years ago

change ContourColor? default back to Paired even if it is 12 color discrete in new matplotlib
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
135C            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          dDELTAdA = (1. - xp(k+it))*dEmindA ! N. B. symm w/r A,B
206          sinDELTA = sind(Delta)
207          cosDELTA = cosd(Delta)
208          RcosDELTA = 1. / cosDELTA
209          tanDELTA = tand(Delta)
210          cosDELTA2 = cosDELTA*cosDELTA     
211          if ( ttheta.le.4500.0 .or. ttheta.gt.13500.0 ) then
212            tmp = sin2THETA2-sinDelta**2
213          else
214            tmp = cosDELTA2-cos2THETA2
215          end if
216          if (tmp .gt. 0) then
217            tmp1 = 1.0/SQRT(tmp)
218            F = abs(cos2THETA) * tmp1
219            dFdA = cosDELTA*cos2THETA*sinDELTA*dDELTAdA*tmp1**3
220          else
221            F = 1.0
222            dFdA = 0.0
223          endif
224c       
225c Calculate G(Delta,2theta) [G = W /(h cos(delta) ] [ FCJ eq. 7(a) and 7(b) ]
226c
227          if ( ttheta.gt.9000.0 ) then                  !are you kidding!! this worked
228!          if(abs(delta-emin) .gt. abs(einfl-emin))then
229              G = 2.0*A*F*RcosDELTA
230              dGdA = 2.0*A*RcosDELTA*(dFdA + F*tanDELTA*dDELTAdA)
231          else
232            G = (-1.0 + ApB*F) * RcosDELTA
233            dGdA = RcosDELTA*(F - tanDELTA*dDELTAdA
234     1             + ApB*F*tanDELTA*dDELTAdA + ApB*dFdA)
235          endif
236
237          WG = wp(k+it) * G
238          sumWG = sumWG + WG
239          sumWRG = sumWRG + WG * R
240          sumWdGdA = sumWdGdA + wp(k+it) * dGdA
241          sumWRdGdA = sumWRdGdA + wp(k+it) * R * dGdA
242          sumWGdRd2t = sumWGdRd2t + WG * dRdT ! N.B. 1/centidegrees
243          sumWGdRdsig = sumWGdRdsig + WG * dRdS
244          sumWGdRdgam = sumWGdRdgam + WG * dRdG
245          sumWGdRdA = sumWGdRdA + WG * dRdT * dDELTAdA ! N. B. symm w/r A,B
246        enddo
247        RsumWG = 1.0/sumWG
248        RsumWG2 = RsumWG * RsumWG
249        PRFUNC = sumWRG * RsumWG
250
251        dydA = (-(sumWRG*sumWdGdA) +
252     1    sumWG*(sumWRdGdA- 100.0 * todeg * sumWGdRdA)) * RsumWG2
253        sigpart = sumWGdRdsig * RsumWG
254        gampart = sumWGdRdgam * RsumWG
255        DPRDT = -SumWGdRd2T * RsumWG
256      else
257c
258C no asymmetry -- nice and simple!
259c
260        CALL PSVOIGT(DTT,SIG,GAM,R,dRdT,dRdS,dRdG)
261        PRFUNC = R
262        dydA = 0.004 * sign(1.0,TTHETA - DTT)
263        sigpart = dRdS
264        gampart = dRdG
265        DPRDT = -dRdT
266      END IF
267      SHLPART = DYDA
268     
269      RETURN
270      END
271     
Note: See TracBrowser for help on using the repository browser.