source:sasview/src/sas/models/c_extension/libigor/2Y_PairCorrelation.c@79492222

Last change on this file since 79492222 was 79492222, checked in by krzywon, 9 years ago

Changed the file and folder names to remove all SANS references.

• Property mode set to `100644`
File size: 4.8 KB
Line
1/*
2 *  PairCorrelation.c
3 *  twoyukawa
4 *
5 *  Created by Marcus Hennig on 5/9/10.
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,error,k;
115        double alpha,real,imag;
116        double Pi = 3.14159265358979323846264338327950288;   /* pi */
117
118
119        for ( n = 0; n < N; n++ ) {
120                data[2*n] = n * ( Sq[n] - 1 );
121                data[2*n+1] = 0;
122        }
123        //      printf("start of new fft\n");
124
125        // data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 })
126        //      int error  = gsl_fft_real_radix2_transform( data, stride, N );
127        error  = 1;
128        dfour1( data-1, N, 1 );         //N is the number of complex points
129
130        //      printf("dfour1 is done\n");
131
132        // if no errors detected
133        if ( error == 1 )
134        {
135                alpha = N * pow( dq, 3 ) / ( 24 * Pi * Pi * phi );
136
137                *dr = 2 * Pi / ( N * dq );
138                for ( k = 0; k < N; k++ )
139                {
140                        // the solutions of the transform is stored in data,
141                        // consult GSL manual for more details
142                        if ( 2*k == 0 || 2*k == 2*N / 2)
143                        {
144                                real = data[2*k];
145                                imag = 0;
146                        }
147                        else if ( 2*k < 2*N / 2 )
148                        {
149                                real = data[2*k];
150                                imag = data[2*k+1];
151                        }
152                        else if ( 2*k > 2*N / 2 )
153                        {
154                                real =  data[2*k];
155                                imag = -data[2*k+1];
156                        }
157
158                        if ( k == 0 )
159                                gr[k] = 0;
160                        else
161//                              gr[k] = 1. + alpha / k * (-imag);               //if using GSL
162                                gr[k] = 1. + alpha / k * (imag);                //if using NR
163                }
164        }
165
166        // release allocated memory
167        free( data );
168
169//      printf(" done with FFT assignment -- Using Numerical Recipes, not GSL\n");
170
171        // return error value
172        return error;
173}
174
175
176// isign == 1 means no scaling of output. isign == -1 multiplies output by nn
177//
178//
179#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
180
181void dfour1(double data[], unsigned long nn, int isign)
182{
183        unsigned long n,mmax,m,j,istep,i;
184        double wtemp,wr,wpr,wpi,wi,theta;
185        double tempr,tempi;
186
187        n=nn << 1;
188        j=1;
189        for (i=1;i<n;i+=2) {
190                if (j > i) {
191                        SWAP(data[j],data[i]);
192                        SWAP(data[j+1],data[i+1]);
193                }
194                m=n >> 1;
195                while (m >= 2 && j > m) {
196                        j -= m;
197                        m >>= 1;
198                }
199                j += m;
200        }
201        mmax=2;
202        while (n > mmax) {
203                istep=mmax << 1;
204                theta=isign*(6.28318530717959/mmax);
205                wtemp=sin(0.5*theta);
206                wpr = -2.0*wtemp*wtemp;
207                wpi=sin(theta);
208                wr=1.0;
209                wi=0.0;
210                for (m=1;m<mmax;m+=2) {
211                        for (i=m;i<=n;i+=istep) {
212                                j=i+mmax;
213                                tempr=wr*data[j]-wi*data[j+1];
214                                tempi=wr*data[j+1]+wi*data[j];
215                                data[j]=data[i]-tempr;
216                                data[j+1]=data[i+1]-tempi;
217                                data[i] += tempr;
218                                data[i+1] += tempi;
219                        }
220                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
221                        wi=wi*wpr+wtemp*wpi+wi;
222                }
223                mmax=istep;
224        }
225}
226#undef SWAP
Note: See TracBrowser for help on using the repository browser.