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 | /* |
---|
45 | int 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 | // |
---|
111 | int 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 | |
---|
179 | void 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 |
---|