| 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 |
|---|
| 73 | c |
|---|
| 74 | c Local Variables for Gaussian Integration |
|---|
| 75 | c |
|---|
| 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/ |
|---|
| 102 | C 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 | |
|---|
| 110 | c |
|---|
| 111 | C f(2theta) intermediates |
|---|
| 112 | c |
|---|
| 113 | TTHETAD = 0.01 * TTHETA |
|---|
| 114 | SIN2THETA = SIND(TTHETAD) |
|---|
| 115 | cos2THETA = COSD(TTHETAD) |
|---|
| 116 | sin2theta2 = sin2THETA * sin2THETA |
|---|
| 117 | cos2THETA2 = cos2THETA * cos2THETA |
|---|
| 118 | c |
|---|
| 119 | C Asymmetry terms |
|---|
| 120 | c |
|---|
| 121 | A = SHL/2. ! A = S/L in FCJ |
|---|
| 122 | ApB = SHL |
|---|
| 123 | ApB2 = ApB*ApB |
|---|
| 124 | c |
|---|
| 125 | C handle the case where there is asymmetry |
|---|
| 126 | c |
|---|
| 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 | |
|---|
| 149 | c |
|---|
| 150 | c Compute Gaussian Quadrature interval |
|---|
| 151 | C Note that the number of points must be large enough so that the |
|---|
| 152 | c interval between 2phi(min) and 2theta must be divided into steps no |
|---|
| 153 | c larger than 0.005 degrees. LWF 8/10/95 |
|---|
| 154 | c |
|---|
| 155 | c Determine which Gauss-Legendre Table to use |
|---|
| 156 | c |
|---|
| 157 | arraynum = 1 |
|---|
| 158 | c |
|---|
| 159 | c Calculate desired number of intervals |
|---|
| 160 | c |
|---|
| 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 |
|---|
| 167 | c |
|---|
| 168 | C Find the next largest set |
|---|
| 169 | c |
|---|
| 170 | do while (arraynum .lt. numtab .and. k.gt.nterms(arraynum)) |
|---|
| 171 | arraynum = arraynum + 1 |
|---|
| 172 | enddo |
|---|
| 173 | ngt = nterms(arraynum) |
|---|
| 174 | c |
|---|
| 175 | C calculate the terms, if they have not been used before |
|---|
| 176 | c |
|---|
| 177 | if (.not. calclfg(arraynum)) then |
|---|
| 178 | calclfg(arraynum) = .true. |
|---|
| 179 | call gauleg(-1.,1.,xpt,wpt,ngt) |
|---|
| 180 | it = fstterm(arraynum)-ngt/2 |
|---|
| 181 | c |
|---|
| 182 | C copy the ngt/2 terms from our working array to the stored array |
|---|
| 183 | c |
|---|
| 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. |
|---|
| 198 | c |
|---|
| 199 | c Compute Convolution integral for 2phi(min) <= delta <= 2theta |
|---|
| 200 | c |
|---|
| 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 |
|---|
| 225 | c |
|---|
| 226 | c Calculate G(Delta,2theta) [G = W /(h cos(delta) ] [ FCJ eq. 7(a) and 7(b) ] |
|---|
| 227 | c |
|---|
| 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 |
|---|
| 258 | c |
|---|
| 259 | C no asymmetry -- nice and simple! |
|---|
| 260 | c |
|---|
| 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 | |
|---|
