source: trunk/fsource/powsubs/epsvoigt.for @ 1484

Last change on this file since 1484 was 1484, checked in by vondreele, 7 years ago

change lower limit of LGmix to 0 from 0.1 - recompile Win-32 & Win-64 fortran sources
fix ellipseSizeDeriv - didn't do it correctly before
fix mustrain & size fxns & derivatives for CW (minor) & TOF major; all checked against numerical derivatives for both lorentzian & gaussian contributions as well as mixing coeff.

File size: 4.3 KB
Line 
1      SUBROUTINE EPSVOIGT(DT,ALP,BET,SIG,GAM,FUNC,DFDX,DFDA,DFDB,
2     1  DFDS,DFDG)
3
4!PURPOSE: Compute function & derivatives exponential X pseudovoigt
5!P.Tompson, D.E. Cox & J.B. Hastings (1987). J. Appl. Cryst.,20,79-83
6
7      INCLUDE       '../INCLDS/COPYRIGT.FOR' 
8
9!PSEUDOCODE:
10
11!CALLING ARGUMENTS:
12
13      REAL*4        DT                  !Delta-TOF from center
14      REAL*4        ALP                 !Exponential rise
15      REAL*4        BET                 !Exponential decay
16      REAL*4        SIG                 !Gaussian variance
17      REAL*4        GAM                 !Lorentzian FWHM
18      REAL*4        FUNC                !Value of pseudo-Voigt at DX
19      REAL*4        DFDX                !dF/dta
20      REAL*4        DFDA                !dF/da
21      REAL*4        DFDB                !dF/db
22      REAL*4        DFDS                !dF/ds
23      REAL*4        DFDG                !dF/dg
24
25!INCLUDE STATEMENTS:
26
27!LOCAL VARIABLES:
28
29      REAL*4        COFG(6),COFL(6)     !Linear combination coeffs
30      REAL*4        ACOFG(7),ACOFL(7)   
31      REAL*4        GNORM               !Gaussian Normalization constant
32      REAL*4        COFT(6),COFN(3)     
33      REAL*4        NORM                !Normalization constant
34      REAL*4        RXA,IXA,RXB,IXB     !Exp-integral arguements
35      REAL*4        RFA,IFA,RFB,IFB     !Exp-integral results
36
37!SUBROUTINES CALLED:
38
39!FUNCTION DEFINITIONS:
40
41!DATA STATEMENTS:
42
43      DATA STOFW/2.35482005/            !2*SQRT(2LN2)
44      DATA SQ2PI/2.506628275/            !SQRT(2PI)
45      DATA PI/3.141592654/            !PI
46      DATA COFT/1.0,2.69269,2.42843,4.47163,0.07842,1.0/
47      DATA COFN/1.36603,-0.47719,0.11116/
48
49!CODE:
50
51      SQSGT = MAX(SQRT(SIG),0.00001)
52      GAM = MAX(GAM,0.00001)
53      FWHG = STOFW*SQSGT
54      PGL = FWHG**5
55      SUMHM = PGL
56      DSDL = 0.0
57      DSDG = 0.0
58      DO ITRM=1,5
59        PGL = PGL/FWHG
60        DSDL = DSDL+FLOAT(ITRM)*COFT(ITRM+1)*PGL
61        DSDG = DSDG+FLOAT(6-ITRM)*COFT(ITRM)*PGL
62        PGL = PGL*GAM
63        SUMHM = SUMHM+COFT(ITRM+1)*PGL
64      END DO
65      FWHM = EXP(0.2*LOG(SUMHM))
66      FRAC = GAM/FWHM
67      DEDF = 0.0
68      PF = 1.0
69      ETA = 0.0
70      DO ITRM=1,3
71        DEDF = DEDF+FLOAT(ITRM)*COFN(ITRM)*PF
72        PF = PF*FRAC
73        ETA = ETA+COFN(ITRM)*PF
74      END DO
75
76      NORM = 0.5*ALP*BET/(ALP+BET)
77      SIGP = (FWHM/STOFW)**2
78      TJ1 = ALP*SIGP
79      TJ2 = BET*SIGP
80      SQSG = SQRT(2.0*SIGP)
81      S2SG = 2.0*SQSG*SIGP
82      Y1 = (TJ1+DT)/SQSG
83      DY1DS = (TJ1-DT)/S2SG
84      Y2 = (TJ2-DT)/SQSG
85      DY2DS = (TJ2+DT)/S2SG
86      EX1 = -0.5*(DT**2)/SIGP
87      DEXDS = -EX1/SIGP
88      H1 = HFUNC(EX1,Y1,0,DH1DY)
89      H2 = HFUNC(EX1,Y2,0,DH2DY)
90
91      TG = NORM*(H1+H2)
92
93      DTGDA = 2.0*NORM*TG/(ALP*ALP)
94      DTGDA = DTGDA+0.5*NORM*SQSG*DH1DY                  !OK
95      DTGDB = 2.0*NORM*TG/(BET*BET)
96      DTGDB = DTGDB+0.5*NORM*SQSG*DH2DY                  !OK
97      DTGDFW = TG*DEXDS+NORM*(DH1DY*DY1DS+DH2DY*DY2DS)
98      DTGDFW = 2.0*DTGDFW*FWHM/STOFW**2                  !OK
99      DTGDT = -NORM*(ALP*H1-BET*H2)                        !OK
100
101      RXA = -ALP*DT
102      RXB = -BET*DT
103      IXA = ALP*FWHM/2.0
104      IXB = BET*FWHM/2.0
105      CALL EXPINT(RXA,IXA,RFA,IFA)
106      CALL EXPINT(RXB,IXB,RFB,IFB)
107      TL = -2.0*NORM*(IFA+IFB)/PI
108      DIVSOR = DT**2+FWHM**2/4.0
109      DTLDT = -2.0*NORM*(ALP*IFA+BET*IFB+FWHM/DIVSOR)/PI      !OK
110
111      DTLDFW = NORM*(-ALP*RFA-BET*RFB-2.0*DT/DIVSOR)/PI
112
113      DTLDA = 2.0*NORM*TL/ALP**2
114      DTLDA = DTLDA+2.0*NORM*(DT*IFA-FWHM*RFA/2.0)/PI            !OK
115      DTLDB = 2.0*NORM*TL/BET**2
116      DTLDB = DTLDB+2.0*NORM*(DT*IFB-FWHM*RFB/2.0)/PI            !OK
117
118      FUNC = ETA*TL+(1.0-ETA)*TG                        !OK
119
120      DFDX = ETA*DTLDT+(1.0-ETA)*DTGDT
121
122      DFWDG = 0.2*DSDG*FWHM/SUMHM
123      DFRDG = -FRAC*DFWDG/FWHM
124
125      DFDS = DEDF*DFRDG*(TL-TG)
126      DFDS = DFDS+(ETA*DTLDFW+(1.0-ETA)*DTGDFW)*DFWDG
127      DFDS = 0.5*DFDS*STOFW/SQSGT                        !OK?
128
129      DFWDL = 0.2*DSDL*FWHM/SUMHM
130      DFRDL = (1.0-FRAC*DFWDL)/FWHM
131
132      DFDG = DEDF*DFRDL*(TL-TG)
133      DFDG = DFDG+(ETA*DTLDFW+(1.0-ETA)*DTGDFW)*DFWDL            !WRONG
134
135      DFDA = ETA*DTLDA+(1.0-ETA)*DTGDA                  !OK
136
137      DFDB = ETA*DTLDB+(1.0-ETA)*DTGDB                  !OK
138
139      RETURN
140      END
Note: See TracBrowser for help on using the repository browser.