source: sans/XOP_Dev/SANSAnalysis/lib/2Y_OneYukawa.c @ 756

Last change on this file since 756 was 756, checked in by srkline, 12 years ago

Adding Yun Liu's 2-Yukawa structure factor to both the library and the XOP. Ideally it would be consolidated to the libStructureFactor files, but there are a lot of files, and labeling them as such seems less confusing.

Currently functions correctly on Mac, still needs to be compiled on Win.

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