source: sasview/sansmodels/src/libigor/2Y_OneYukawa.c @ b11e127

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since b11e127 was a42ea15, checked in by Jae Cho <jhjcho@…>, 13 years ago

Yukawa structure model kernels (IGOR/NIST)
  • Property mode set to 100644
File size: 10.6 KB
RevLine 
[a42ea15]1/*
2 *  Yukawa.c
3 *  twoyukawa
4 *
5 *  Created by Marcus Hennig on 5/12/10.
6 *  Copyright 2010 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10#include "2Y_OneYukawa.h"
11#include "2Y_cpoly.h"
12#include "2Y_utility.h"
13#include "2Y_PairCorrelation.h"
14#include <stdio.h>
15#include <math.h>
16#include <stdlib.h>
17
18/*
19==================================================================================================
20 
21 Structure factor for the potential
22 
23 V(r) = -kB * T * K * exp[ -Z * (r - 1)] / r for r > 1
24 V(r) = inf for r <= 1
25 
26 The structure factor is parametrized by (a, b, c, d)
27 which depend on (Z, K, phi). 
28 
29==================================================================================================
30*/
31
32double Y_sigma( double s, double Z, double a, double b, double c, double d )
33{
34        return -(a / 2. + b + c * exp( -Z )) / s + a * pow( s, -3 ) + b * pow( s, -2 ) + ( c + d ) * pow( s + Z, -1 );
35}
36
37double Y_tau( double s, double Z, double a, double b, double c )
38{
39        return b * pow( s, -2 ) + a * ( pow( s, -3 ) + pow( s, -2 ) ) - pow( s, -1 ) * c * Z * exp( -Z ) * pow( s + Z, -1 );
40} 
41
42double Y_q( double s, double Z, double a, double b, double c, double d )
43{
44        return Y_sigma(s, Z, a, b, c, d ) - exp( -s ) * Y_tau( s, Z, a, b, c );
45}
46
47double Y_g( double s, double phi, double Z, double a, double b, double c, double d )
48{
49        return s * Y_tau( s, Z, a, b, c ) * exp( -s ) / ( 1 - 12 * phi * Y_q( s, Z, a, b, c, d ) );
50}
51
52double Y_hq( double q, double Z, double K, double v )
53{
54        double t1, t2, t3, t4;
55
56        if ( q == 0) 
57        {
58                return (exp(-2*Z)*(v + (v*(-1 + Z) - 2*K*Z)*exp(Z))*(-(v*(1 + Z)) + (v + 2*K*Z*(1 + Z))*exp(Z))*pow(K,-1)*pow(Z,-4))/4.;
59        }
60        else 
61        {
62               
63                t1 = ( 1 - v / ( 2 * K * Z * exp( Z ) ) ) * ( ( 1 - cos( q ) ) / ( q*q ) - 1 / ( Z*Z + q*q ) );
64                t2 = ( v*v * ( q * cos( q ) - Z * sin( q ) ) ) / ( 4 * K * Z*Z * q * ( Z*Z + q*q ) );
65                t3 = ( q * cos( q ) + Z * sin( q ) ) / ( q * ( Z*Z + q*q ) );
66                t4 = v / ( Z * exp( Z ) ) - v*v / ( 4 * K * Z*Z * exp( 2 * Z ) ) - K;
67               
68                return v / Z * t1 - t2 + t3 * t4;
69        }
70}
71
72double Y_pc( double q, 
73                        double Z, double K, double phi,
74                        double a, double b, double c, double d )
75{       
76        double v = 24 * phi * K * exp( Z ) * Y_g( Z, phi, Z, a, b, c, d );
77       
78        double a0 = a * a;
79        double b0 = -12 * phi *( pow( a + b,2 ) / 2 + a * c * exp( -Z ) );
80       
81        double t1, t2, t3, t4;
82       
83        if ( q == 0 ) 
84        {
85                t1 = a0 / 3;
86                t2 = b0 / 4;
87                t3 = a0 * phi / 12;
88        }
89        else 
90        {
91                t1 = a0 * ( sin( q ) - q * cos( q ) ) / pow( q, 3 );
92                t2 = b0 * ( 2 * q * sin( q ) - ( q * q - 2 ) * cos( q ) - 2 ) / pow( q, 4 );
93                t3 = a0 * phi * ( ( q*q - 6 ) * 4 * q * sin( q ) - ( pow( q, 4 ) - 12 * q*q + 24) * cos( q ) + 24 ) / ( 2 * pow( q, 6 ) );
94        }
95        t4 = Y_hq( q, Z, K, v );
96        return -24 * phi * ( t1 + t2 + t3 + t4 );
97}
98
99double SqOneYukawa( double q, 
100                                 double Z, double K, double phi,
101                                 double a, double b, double c, double d )
102{
103        //structure factor one-yukawa potential
104        return 1. / ( 1. - Y_pc( q, Z, K, phi, a, b, c, d ) );
105}
106
107
108/*
109==================================================================================================
110
111 The structure factor S(q) for the one-Yukawa potential is parameterized by a,b,c,d which are
112 the solution of a system of non-linear equations given Z,K, phi. There are at most 4 solutions
113 from which the physical one is chosen
114 
115==================================================================================================
116 */
117
118double Y_LinearEquation_1( double Z, double K, double phi, double a, double b, double c, double d )
119{
120        return b - 12*phi*(-a/8. - b/6. + d*pow(Z,-2) + c*(pow(Z,-2) - exp(-Z)*(0.5 + (1 + Z)*pow(Z,-2))));
121}
122
123double Y_LinearEquation_2( double Z, double K, double phi, double a, double b, double c, double d )
124{
125        return 1 - a - 12*phi*(-a/3. - b/2. + d*pow(Z,-1) + c*(pow(Z,-1) - (1 + Z)*exp(-Z)*pow(Z,-1)));
126}
127
128double Y_LinearEquation_3( double Z, double K, double phi, double a, double b, double c, double d )
129{
130        return K*exp(Z) - Z*d*(1-12*phi*Y_q(Z, Z, a, b, c, d));
131}
132
133double Y_NonlinearEquation( double Z, double K, double phi, double a, double b, double c, double d )
134{
135        return c + d - 12*phi*((c + d)*Y_sigma(Z, Z, a, b, c, d) - c*exp(-Z)*Y_tau(Z, Z, a, b, c));
136}
137
138// Check the computed solutions satisfy the system of equations
139int Y_CheckSolution( double Z, double K, double phi, 
140                                         double a, double b, double c, double d )
141{
142        double eq_1 = chop( Y_LinearEquation_1 ( Z, K, phi, a, b, c, d ) );
143        double eq_2 = chop( Y_LinearEquation_2 ( Z, K, phi, a, b, c, d ) );
144        double eq_3 = chop( Y_LinearEquation_3 ( Z, K, phi, a, b, c, d ) );
145        double eq_4 = chop( Y_NonlinearEquation( Z, K,  phi, a, b, c, d ) );
146       
147        // check if all equation are zero
148        return eq_1 == 0 && eq_2 == 0 && eq_3 == 0 && eq_4 == 0;
149}
150
151int Y_SolveEquations( double Z, double K, double phi, double* a, double* b, double* c, double* d, int debug )
152{
153        char buf[256];
154
155        // at most there are 4 solutions for a,b,c,d
156        double sol_a[4], sol_b[4], sol_c[4], sol_d[4];
157       
158        double m11 = (3*phi)/2.;
159        double m13 = 6*phi*exp(-Z)*(2 + Z*(2 + Z) - 2*exp(Z))*pow(Z,-2);
160        double m23 = -12*phi*exp(-Z)*(-1 - Z + exp(Z))*pow(Z,-1);
161        double m31 = -6*phi*exp(-Z)*pow(Z,-2)*(2*(1 + Z) + exp(Z)*(-2 + pow(Z,2)));
162        double m32 = -12*phi*(-1 + Z + exp(-Z))*pow(Z,-1);
163        double m33 = 6*phi*exp(-2*Z)*pow(-1 + exp(Z),2);
164       
165        double delta = m23*m31 - m13*m32 + m11*(-4*m13*m31 + (4*m23*m31)/3. + (8*m13*m32)/3. - m23*m32) + m33 + (4*(-3 + m11)*m11*m33)/9.;
166        double a1 = -(K*(m23 + (4*m11*(-3*m13 + m23))/3.)*exp(Z));
167        double a2 = -(m13*(m32 + 4*m11*Z)) + ((3 + 4*m11)*(m33 + m23*Z))/3.;
168        double a3 = -2*phi*pow(Z,-2)*(6*m23*m32 - 24*m11*m33 + 2*Z*((3 + 4*m11)*m33 - 3*m13*(m32 + 2*m11*Z)) + (3 + 4*m11)*m23*pow(Z,2));
169       
170        double b1 = -(K*((-3 + 8*m11)*m13 - 3*m11*m23)*exp(Z))/3.;
171        double b2 = m13*(m31 - Z + (8*m11*Z)/3.) - m11*(m33 + m23*Z);
172        double b3 = 2*phi*pow(Z,-2)*(m13*Z*(-6*m31 + 3*Z - 8*m11*Z) + 2*m33*(3 - 8*m11 + 3*m11*Z) + 3*m23*(2*m31 + m11*pow(Z,2)));
173       
174        double c1 = -(K*exp(Z)*pow(3 - 2*m11,2))/9.;
175        double c2 = -((3 + 4*m11)*m31)/3. + m11*m32 + Z + (4*(-3 + m11)*m11*Z)/9.;
176        double c3 = (-2*phi*pow(Z,-2)*(6*(12*m11*m31 + 3*m32 - 8*m11*m32) - 6*((3 + 4*m11)*m31 - 3*m11*m32)*Z + pow(3 - 2*m11,2)*pow(Z,2)))/3.;
177       
178        // determine d, as roots of the polynomial, from that build a,b,c
179       
180        double real_coefficient[5];
181        double imag_coefficient[5];
182       
183        double real_root[4];
184        double imag_root[4];
185       
186        double zeta = 24*phi*pow(-6*phi*Z*cosh(Z/2.) + (12*phi + (-1 + phi)*pow(Z,2))*sinh(Z/2.),2);   
187        double A[5];
188        int degree,i,j,n_roots;
189        double x,y;
190        int n,selected_root;
191        double qmax,q,dq,min,sum,dr;
192        double *sq,*gr;
193       
194       
195        A[0] = -(exp(3*Z)*pow(K,2)*pow(-1 + phi,2)*pow(Z,3) / zeta );
196        A[1] = K*Z*exp(Z)*(6*phi*(2 + 4*phi + (2 + phi)*Z) + exp(Z)*
197                                           ((-24 + Z*(18 + (-6 + Z)*Z))*pow(phi,2) - 2*phi*(6 + (-3 + Z)*pow(Z,2)) + pow(Z,3))) / zeta;
198        A[2] = -12*K*phi*exp(Z)*pow(-1 + phi,2)*pow(Z,3)/zeta;
199        A[3] = 6*phi*Z*exp(-Z)*(-12*phi*(1 + 2*phi)*(-1 + exp(Z)) + 6*phi*Z*(3*phi + (2 + phi)*exp(Z)) + 
200                                                        6*(-1 + phi)*phi*pow(Z,2) + pow(-1 + phi,2)*pow(Z,3))/zeta;
201        A[4] = -36*exp(-Z)*pow(-1 + phi,2)*pow(phi,2)*pow(Z,3)/zeta;
202/*     
203        if ( debug )
204        {
205                sprintf (buf, "(Z,K,phi) = (%g, %g, %g)\r", K, Z, phi );
206                XOPNotice(buf);
207                sprintf (buf, "A = (%g, %g, %g, %g, %g)\r", A[0], A[1], A[2], A[3], A[4] );
208                XOPNotice(buf);
209        }
210*/     
211        //integer degree of polynomial
212        degree = 4;
213       
214        // vector of real and imaginary coefficients in order of decreasing powers
215        for ( i = 0; i <= degree; i++ )
216        {
217                real_coefficient[i] = A[4-i];
218                imag_coefficient[i] = 0.;
219        }
220       
221        // Jenkins-Traub method to approximate the roots of the polynomial
222        cpoly( real_coefficient, imag_coefficient, degree, real_root, imag_root );
223       
224        // show the result if in debug mode
225/*
226        if ( debug )
227        {
228                for ( i = 0; i < degree; i++ )
229                {
230                        x = real_root[i];
231                        y = imag_root[i];
232                        sprintf(buf, "root(%d) = %g + %g i\r", i+1, x, y );
233                        XOPNotice(buf);
234                }
235                sprintf(buf, "\r" );
236                XOPNotice(buf);
237        }
238 */
239        // determine the set of solutions for a,b,c,d,
240        j = 0;
241        for ( i = 0; i < degree; i++ ) 
242        {
243                x = real_root[i];
244                y = imag_root[i];
245               
246                if ( chop( y ) == 0 ) 
247                {
248                        sol_a[j] = ( a1 + a2 * x + a3 * x * x ) / ( delta * x );
249                        sol_b[j] = ( b1 + b2 * x + b3 * x * x ) / ( delta * x );
250                        sol_c[j] = ( c1 + c2 * x + c3 * x * x ) / ( delta * x );
251                        sol_d[j] = x;
252       
253                        j++;
254                }
255        }
256
257        // number  remaining roots
258        n_roots = j;
259       
260//      sprintf(buf, "inside Y_solveEquations OK, before malloc: n_roots = %d\r",n_roots);
261//      XOPNotice(buf);
262       
263        // if there is still more than one root left, than choose the one with the minimum
264        // average value inside the hardcore
265       
266       
267       
268        if ( n_roots > 1 )
269        {
270                // the number of q values should be a power of 2
271                // in order to speed up the FFT
272                n = 1 << 14;            //= 16384
273               
274                // the maximum q value should be large enough
275                // to enable a reasoble approximation of g(r)
276                qmax = 1000.;
277                dq = qmax / ( n - 1 );
278               
279                // step size for g(r) = dr
280               
281                // allocate memory for pair correlation function g(r)
282                // and structure factor S(q)
283                // (note that sq and gr are pointers)
284                sq = malloc( sizeof( double ) * n );
285                gr = malloc( sizeof( double ) * n );
286               
287                // loop over all remaining roots
288                min = 1e99;
289                selected_root=0;       
290               
291//              sprintf(buf, "after malloc: n,dq = %d  %g\r",n,dq);
292//              XOPNotice(buf);
293                               
294                for ( j = 0; j < n_roots; j++) 
295                {
296                        // calculate structure factor at different q values
297                        for ( i = 0; i < n ; i++) 
298                        {
299                                q = dq * i;
300                                sq[i] = SqOneYukawa( q, Z, K, phi, sol_a[j], sol_b[j], sol_c[j], sol_d[j] );
301/*                             
302                                if(i<20 && debug) {
303                                        sprintf(buf, "after SqOneYukawa: s(q) = %g\r",sq[i] );
304                                        XOPNotice(buf);
305                                }
306 */
307                        }
308                                               
309                        // calculate pair correlation function for given
310                        // structure factor, g(r) is computed at values
311                        // r(i) = i * dr
312                        PairCorrelation( phi, dq, sq, &dr, gr, n );
313                       
314//                      sprintf(buf, "after PairCorrelation: \r");
315//                      XOPNotice(buf);
316                       
317                        // determine sum inside the hardcore
318                        // 0 =< r < 1 of the pair-correlation function
319                        sum = 0;
320                        for (i = 0; i < floor( 1. / dr ); i++ ) 
321                        {
322                                sum += fabs( gr[i] );
323/*                             
324                                if(i<20 && debug) {
325                                        sprintf(buf, "g(r) in core = %g\r",fabs(gr[i]));
326                                        XOPNotice(buf);
327                                }
328 */
329                        }
330                       
331//                      sprintf(buf, "after hard core: sum, min = %g %g\r",sum,min);
332//                      XOPNotice(buf);
333                       
334                        if ( sum < min )
335                        {
336                                min = sum;
337                                selected_root = j;
338                        }
339                       
340                       
341                }       
342                free( gr );
343                free( sq );
344               
345//              sprintf(buf, "after free: selected root = %d\r",selected_root);
346//              XOPNotice(buf);
347               
348                // physical solution was found
349                *a = sol_a[ selected_root ];
350                *b = sol_b[ selected_root ];
351                *c = sol_c[ selected_root ];
352                *d = sol_d[ selected_root ];
353               
354//              sprintf(buf, "after solution found (ret 1): \r");
355//              XOPNotice(buf);
356               
357                return 1;
358        }
359        else if ( n_roots == 1 ) 
360        {
361                *a = sol_a[0];
362                *b = sol_b[0];
363                *c = sol_c[0];
364                *d = sol_d[0];
365               
366                return 1;
367        }
368        // no solution was found
369        return 0;
370}
Note: See TracBrowser for help on using the repository browser.