source: sans/XOP_Dev/SANSAnalysis/lib/2Y_PairCorrelation.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: 4.8 KB
Line 
1/*
2 *  PairCorrelation.c
3 *  twoyukawa
4 *
5 *  Created by Marcus Hennig on 5/9/10.
6 *  Copyright 2010 __MyCompanyName__. All rights reserved.
7 *
8 */
9
10#include "2Y_PairCorrelation.h"
11#include <stdio.h>
12#include <stdlib.h>
13#include <math.h>
14//#include <gsl/gsl_errno.h>
15//#include <gsl/gsl_fft_real.h>
16
17/*
18===================================================================================================
19
20 Source: J.B.Hayter: A Program for the fast bi-directional transforms
21 between g(r) and S(Q), ILL internal scientific report, October 1979
22 
23 The transformation between structure factor and pair correlation
24 function is given by
25 
26 g(x) = 1 + 1 / (12*pi*phi*x) * int( [S(q)-q]*q*sin(q*x), { 0, inf } )
27 
28 where phi is the volume frcation, x and q are dimensionless variables,
29 scaled by the radius a of the particles:
30 
31 r = x * a;
32 Q = q / a
33 
34 Discretizing the integral leads to
35 
36 x[k] = 2*pi*k / (N * dq)
37 g(x[k]) = 1 + N*dq^3 / (24*pi^2*phi*k) * Im{ sum(S[n]*exp(2*pi*i*n*k/N),{n,0,N-1}) }
38 
39 where S[n] = n*(S(q[n])-1) with q[n]=n * dq
40
41===================================================================================================
42*/
43
44/*
45int PairCorrelation_GSL( double phi, double dq, double* Sq, double* dr, double* gr, int N )
46{
47        double* data = malloc( sizeof(double) * N );
48        int n;
49        for ( n = 0; n < N; n++ )
50                data[n] = n * ( Sq[n] - 1 );
51       
52        // data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
53        int stride = 1;
54        int error  = gsl_fft_real_radix2_transform( data, stride, N );
55       
56        // if no errors detected
57        if ( error == GSL_SUCCESS )
58        {
59                double alpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi );
60       
61       
62                *dr = 2 * M_PI / ( N * dq ); 
63                int k;
64                double real, imag;
65                for ( k = 0; k < N; k++ )
66                {
67                        // the solutions of the transform is stored in data,
68                        // consult GSL manual for more details
69                        if ( k == 0 || k == N / 2)
70                        {
71                                real = data[k];
72                                imag = 0;
73                        }
74                        else if ( k < N / 2 )
75                        {
76                                real = data[k];
77                                imag = data[N-k];
78                        }
79                        else if ( k > N / 2 )
80                        {
81                                real =  data[N-k];
82                                imag = -data[k];       
83                        }
84                       
85                        if ( k == 0 )
86                                gr[k] = 0;
87                        else
88                                gr[k] = 1. + alpha / k * (-imag);
89                }
90        }
91        // if N is not a power of two
92        else if ( error == GSL_EDOM )
93        {
94                printf( "N is not a power of 2\n" );
95        }
96        else
97        {
98                printf( "Could not perform DFT (discrete fourier transform)\n" );
99        }
100        // release allocated memory
101        free( data );
102       
103        // return error value
104        return error;
105}
106 */
107
108
109// this uses numerical recipes for the FFT
110//
111int PairCorrelation( double phi, double dq, double* Sq, double* dr, double* gr, int N )
112{
113        double* data = malloc( sizeof(double) * N * 2);
114        int n;
115        for ( n = 0; n < N; n++ ) {
116                data[2*n] = n * ( Sq[n] - 1 );
117                data[2*n+1] = 0;
118        }
119        //      printf("start of new fft\n");
120       
121        // data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
122        //      int error  = gsl_fft_real_radix2_transform( data, stride, N );
123        int error  = 1;
124        dfour1( data-1, N, 1 );         //N is the number of complex points
125       
126        //      printf("dfour1 is done\n");
127       
128        // if no errors detected
129        if ( error == 1 ) 
130        {
131                double alpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi );
132               
133                *dr = 2 * M_PI / ( N * dq ); 
134                int k;
135                double real, imag;
136                for ( k = 0; k < N; k++ )
137                {
138                        // the solutions of the transform is stored in data,
139                        // consult GSL manual for more details
140                        if ( 2*k == 0 || 2*k == 2*N / 2)
141                        { 
142                                real = data[2*k];
143                                imag = 0;
144                        }
145                        else if ( 2*k < 2*N / 2 ) 
146                        {
147                                real = data[2*k];
148                                imag = data[2*k+1];
149                        }
150                        else if ( 2*k > 2*N / 2 )
151                        {
152                                real =  data[2*k];
153                                imag = -data[2*k+1];   
154                        }
155                       
156                        if ( k == 0 ) 
157                                gr[k] = 0;
158                        else
159//                              gr[k] = 1. + alpha / k * (-imag);               //if using GSL
160                                gr[k] = 1. + alpha / k * (imag);                //if using NR
161                } 
162        } 
163       
164        // release allocated memory
165        free( data );
166       
167//      printf(" done with FFT assignment -- Using Numerical Recipes, not GSL\n");
168       
169        // return error value
170        return error;
171}
172
173
174// isign == 1 means no scaling of output. isign == -1 multiplies output by nn
175//
176//
177#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
178
179void dfour1(double data[], unsigned long nn, int isign)
180{
181        unsigned long n,mmax,m,j,istep,i;
182        double wtemp,wr,wpr,wpi,wi,theta;
183        double tempr,tempi;
184       
185        n=nn << 1;
186        j=1;
187        for (i=1;i<n;i+=2) {
188                if (j > i) {
189                        SWAP(data[j],data[i]);
190                        SWAP(data[j+1],data[i+1]);
191                }
192                m=n >> 1;
193                while (m >= 2 && j > m) {
194                        j -= m;
195                        m >>= 1;
196                }
197                j += m;
198        }
199        mmax=2;
200        while (n > mmax) {
201                istep=mmax << 1;
202                theta=isign*(6.28318530717959/mmax);
203                wtemp=sin(0.5*theta);
204                wpr = -2.0*wtemp*wtemp;
205                wpi=sin(theta);
206                wr=1.0;
207                wi=0.0;
208                for (m=1;m<mmax;m+=2) {
209                        for (i=m;i<=n;i+=istep) {
210                                j=i+mmax;
211                                tempr=wr*data[j]-wi*data[j+1];
212                                tempi=wr*data[j+1]+wi*data[j];
213                                data[j]=data[i]-tempr;
214                                data[j+1]=data[i+1]-tempi;
215                                data[i] += tempr;
216                                data[i+1] += tempi;
217                        }
218                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
219                        wi=wi*wpr+wtemp*wpi+wi;
220                }
221                mmax=istep;
222        }
223}
224#undef SWAP
Note: See TracBrowser for help on using the repository browser.