Changeset 760 for sans/XOP_Dev/SANSAnalysis
- Timestamp:
- Oct 27, 2010 9:42:29 AM (12 years ago)
- Location:
- sans/XOP_Dev/SANSAnalysis/lib
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
sans/XOP_Dev/SANSAnalysis/lib/2Y_OneYukawa.c
r756 r760 52 52 double Y_hq( double q, double Z, double K, double v ) 53 53 { 54 double t1, t2, t3, t4; 55 54 56 if ( q == 0) 55 57 { … … 58 60 else 59 61 { 60 double t1, t2, t3, t4;61 62 62 63 t1 = ( 1 - v / ( 2 * K * Z * exp( Z ) ) ) * ( ( 1 - cos( q ) ) / ( q*q ) - 1 / ( Z*Z + q*q ) ); … … 78 79 double b0 = -12 * phi *( pow( a + b,2 ) / 2 + a * c * exp( -Z ) ); 79 80 80 double t1, t2, t3 ;81 double t1, t2, t3, t4; 81 82 82 83 if ( q == 0 ) … … 92 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 ) ); 93 94 } 94 doublet4 = Y_hq( q, Z, K, v );95 t4 = Y_hq( q, Z, K, v ); 95 96 return -24 * phi * ( t1 + t2 + t3 + t4 ); 96 97 } … … 185 186 double zeta = 24*phi*pow(-6*phi*Z*cosh(Z/2.) + (12*phi + (-1 + phi)*pow(Z,2))*sinh(Z/2.),2); 186 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 187 194 188 195 A[0] = -(exp(3*Z)*pow(K,2)*pow(-1 + phi,2)*pow(Z,3) / zeta ); … … 193 200 6*(-1 + phi)*phi*pow(Z,2) + pow(-1 + phi,2)*pow(Z,3))/zeta; 194 201 A[4] = -36*exp(-Z)*pow(-1 + phi,2)*pow(phi,2)*pow(Z,3)/zeta; 195 202 /* 196 203 if ( debug ) 197 204 { … … 201 208 XOPNotice(buf); 202 209 } 203 210 */ 204 211 //integer degree of polynomial 205 intdegree = 4;212 degree = 4; 206 213 207 214 // vector of real and imaginary coefficients in order of decreasing powers 208 int i;209 215 for ( i = 0; i <= degree; i++ ) 210 216 { … … 217 223 218 224 // show the result if in debug mode 219 double x, y; 225 /* 220 226 if ( debug ) 221 227 { … … 230 236 XOPNotice(buf); 231 237 } 238 */ 232 239 // determine the set of solutions for a,b,c,d, 233 intj = 0;240 j = 0; 234 241 for ( i = 0; i < degree; i++ ) 235 242 { … … 249 256 250 257 // number remaining roots 251 intn_roots = j;258 n_roots = j; 252 259 253 260 // sprintf(buf, "inside Y_solveEquations OK, before malloc: n_roots = %d\r",n_roots); … … 256 263 // if there is still more than one root left, than choose the one with the minimum 257 264 // average value inside the hardcore 265 266 267 258 268 if ( n_roots > 1 ) 259 269 { 260 270 // the number of q values should be a power of 2 261 271 // in order to speed up the FFT 262 intn = 1 << 14; //= 16384272 n = 1 << 14; //= 16384 263 273 264 274 // the maximum q value should be large enough 265 275 // 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; 276 qmax = 1000.; 277 dq = qmax / ( n - 1 ); 278 279 // step size for g(r) = dr 271 280 272 281 // allocate memory for pair correlation function g(r) 273 282 // and structure factor S(q) 274 double* sq = malloc( sizeof( double ) * n ); 275 double* gr = malloc( sizeof( double ) * n ); 283 // (note that sq and gr are pointers) 284 sq = malloc( sizeof( double ) * n ); 285 gr = malloc( sizeof( double ) * n ); 276 286 277 287 // loop over all remaining roots 278 doublemin = INFINITY;279 intselected_root=0;288 min = INFINITY; 289 selected_root=0; 280 290 281 291 // sprintf(buf, "after malloc: n,dq = %d %g\r",n,dq); … … 289 299 q = dq * i; 290 300 sq[i] = SqOneYukawa( q, Z, K, phi, sol_a[j], sol_b[j], sol_c[j], sol_d[j] ); 291 301 /* 292 302 if(i<20 && debug) { 293 303 sprintf(buf, "after SqOneYukawa: s(q) = %g\r",sq[i] ); 294 304 XOPNotice(buf); 295 305 } 306 */ 296 307 } 297 308 … … 306 317 // determine sum inside the hardcore 307 318 // 0 =< r < 1 of the pair-correlation function 308 doublesum = 0;319 sum = 0; 309 320 for (i = 0; i < floor( 1. / dr ); i++ ) 310 321 { 311 322 sum += fabs( gr[i] ); 312 323 /* 313 324 if(i<20 && debug) { 314 325 sprintf(buf, "g(r) in core = %g\r",fabs(gr[i])); 315 326 XOPNotice(buf); 316 327 } 328 */ 317 329 } 318 330 -
sans/XOP_Dev/SANSAnalysis/lib/2Y_PairCorrelation.c
r756 r760 112 112 { 113 113 double* data = malloc( sizeof(double) * N * 2); 114 int n; 114 int n,error,k; 115 double alpha,real,imag; 116 115 117 for ( n = 0; n < N; n++ ) { 116 118 data[2*n] = n * ( Sq[n] - 1 ); … … 121 123 // data[k] -> sum( data[n] * exp(-2*pi*i*n*k/N), {n, 0, N-1 }) 122 124 // int error = gsl_fft_real_radix2_transform( data, stride, N ); 123 interror = 1;125 error = 1; 124 126 dfour1( data-1, N, 1 ); //N is the number of complex points 125 127 … … 129 131 if ( error == 1 ) 130 132 { 131 doublealpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi );133 alpha = N * pow( dq, 3 ) / ( 24 * M_PI * M_PI * phi ); 132 134 133 135 *dr = 2 * M_PI / ( N * dq ); 134 int k;135 double real, imag;136 136 for ( k = 0; k < N; k++ ) 137 137 { -
sans/XOP_Dev/SANSAnalysis/lib/2Y_TwoYukawa.c
r756 r760 84 84 double TY_hq( double q, double Z, double K, double v ) 85 85 { 86 double t1, t2, t3, t4; 87 86 88 if ( q == 0) 87 89 { … … 90 92 else 91 93 { 92 double t1, t2, t3, t4;93 94 94 95 t1 = ( 1 - v / ( 2 * K * Z * exp( Z ) ) ) * ( ( 1 - cos( q ) ) / ( q*q ) - 1 / ( Z*Z + q*q ) ); … … 111 112 double b0 = -12 * phi *( pow( a + b,2 ) / 2 + a * ( c1 * exp( -Z1 ) + c2 * exp( -Z2 ) ) ); 112 113 113 double t1, t2, t3 ;114 double t1, t2, t3,t4; 114 115 115 116 if ( q == 0 ) … … 125 126 t3 = a0 * phi * ( ( q*q - 6 ) * 4 * q * sin( q ) - ( pow( q, 4 ) - 12 * q*q + 24) * cos( q ) + 24 ) / ( 2 * pow( q, 6 ) ); 126 127 } 127 doublet4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 );128 t4 = TY_hq( q, Z1, K1, v1 ) + TY_hq( q, Z2, K2, v2 ); 128 129 return -24 * phi * ( t1 + t2 + t3 + t4 ); 129 130 } … … 237 238 (m33*(3*m41 + 4*m11*m41 - 3*m11*m42) + (-3*m31 - 4*m11*m31 + 3*m11*m32)*m43) + 3*m23*(-3*m34*m41 - 4*m11*m34*m41 + 238 239 3*m11*m34*m42 + 3*m31*m44 + 4*m11*m31*m44 - 3*m11*m32*m44) - (m34*m43 - m33*m44)*pow(3 - 2*m11,2))/9.; 239 240 /* 240 241 if( debug ) 241 242 { … … 245 246 XOPNotice(buf); 246 247 } 248 */ 247 249 248 250 /* Matrix representation of the determinant of the of the system where row refering to … … 272 274 2*m11*(-3*(m14*m43 - m13*m44)*(Z1 + Z2)*pow(Z1,2) + 2*m34*(m43*(-3 + Z1)*(Z1 + Z2) - 3*m13*Z2*pow(Z1,2)) + 273 275 m33*(-2*m44*(-3 + Z1)*(Z1 + Z2) + 6*m14*Z2*pow(Z1,2)))))*pow(Z1 + Z2,-1); 274 276 /* 275 277 if( debug ) 276 278 { … … 283 285 XOPNotice(buf); 284 286 } 287 */ 285 288 286 289 /* Matrix representation of the determinant of the of the system where row refering to … … 309 312 m11*Z1*(2*m34*m43*(8 - 3*Z1) + 2*m33*m44*(-8 + 3*Z1) + (8*m14*m43 - 3*m24*m43 - 8*m13*m44 + 3*m23*m44)*pow(Z1,2)))* 310 313 pow(Z1 + Z2,-1) + 6*(-(m14*(m23*m31 + m33)) + m13*(m24*m31 + m34))*pow(Z1,3)*pow(Z1 + Z2,-1)); 311 314 /* 312 315 if( debug ) 313 316 { … … 320 323 XOPNotice(buf); 321 324 } 325 */ 322 326 323 327 /* Matrix representation of the determinant of the of the system where row refering to … … 348 352 Z2*(4*(3*m31 - 2*m32)*m44 + Z1*(-4*m31*m44 + 3*m32*m44 - 2*(m14*(-6*m31 + 4*m32 + 3*m41 - 2*m42) + m44)*Z1) + 349 353 m34*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))))*pow(Z1 + Z2,-1))/3.; 350 354 /* 351 355 if( debug ) 352 356 { … … 359 363 XOPNotice(buf); 360 364 } 365 */ 361 366 /* Matrix representation of the determinant of the of the system where row refering to 362 367 the variable c1 is replaced by solution vector */ … … 387 392 Z2*(4*(3*m31 - 2*m32)*m43 + Z1*(-4*m31*m43 + 3*m32*m43 - 2*(m13*(-6*m31 + 4*m32 + 3*m41 - 2*m42) + m43)*Z1) + 388 393 m33*(m42*(8 - 3*Z1) + 4*m41*(-3 + Z1) + 4*pow(Z1,2)))))*pow(Z1 + Z2,-1))/3.; 389 394 /* 390 395 if( debug ) 391 396 { … … 398 403 XOPNotice(buf); 399 404 } 405 */ 400 406 401 407 /* coefficient matrices of nonlinear equation 1 */ … … 489 495 TY_A43 /= norm_A; 490 496 TY_A52 /= norm_A;*/ 491 497 /* 492 498 if( debug ) 493 499 { … … 505 511 XOPNotice(buf); 506 512 } 513 */ 514 507 515 /* coefficient matrices of nonlinear equation 2 */ 508 516 … … 597 605 TY_B32 /= norm_B; 598 606 TY_B34 /= norm_B; */ 599 607 /* 600 608 if( debug ) 601 609 { … … 609 617 XOPNotice(buf); 610 618 } 619 */ 611 620 612 621 /* decrease order of nonlinear equation 1 by means of equation 2 */ … … 631 640 TY_F39 = 2*TY_A52*TY_B24*TY_B25 - TY_A43*TY_B24*TY_B34 - TY_A42*TY_B25*TY_B34; 632 641 TY_F310 = -(TY_A43*TY_B25*TY_B34) + TY_A52*pow(TY_B25,2); 633 642 643 /* 634 644 if( debug ) 635 645 { … … 643 653 XOPNotice(buf); 644 654 } 645 655 */ 656 646 657 TY_G13 = -(TY_B12*TY_F32); 647 658 TY_G14 = -(TY_B12*TY_F33); … … 668 679 TY_G213 = -(TY_B24*TY_F310) - TY_B25*TY_F39; 669 680 TY_G214 = -(TY_B25*TY_F310); 670 681 /* 671 682 if( debug ) 672 683 { … … 678 689 XOPNotice(buf); 679 690 } 680 691 */ 681 692 // coefficients for polynomial 682 693 TY_w[0] = (-(TY_A21*TY_B12) + TY_A12*TY_B21)*(TY_A52*TY_B21 - TY_A41*TY_B32)*pow(TY_B21,2)*pow(TY_B32,3); … … 828 839 829 840 TY_w[22] = (-(TY_A23*TY_B14) + TY_A12*TY_B25)*(TY_A52*TY_B25 - TY_A43*TY_B34)*pow(TY_B25,2)*pow(TY_B34,3); 830 841 /* 831 842 if( debug ) 832 843 { … … 841 852 XOPNotice(buf); 842 853 } 854 */ 843 855 } 844 856 … … 890 902 int debug ) 891 903 { 892 // reduce system to a polynomial from which all solution are extracted893 // by doing that a lot of global background variables are set894 TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, debug );895 904 896 905 // the two coupled non-linear eqautions were reduced to a … … 907 916 //integer degree of polynomial 908 917 int degree = 22; 918 int i; 919 double x,y; 920 double var_a, var_b, var_c1, var_c2, var_d1, var_d2; 921 double sol_a[22], sol_b[22], sol_c1[22], sol_c2[22], sol_d1[22], sol_d2[22]; 922 923 int j = 0; 924 925 int n_roots,n,selected_root; 926 double dr,qmax,q,dq,min,sum; 927 double *sq,*gr; 928 929 930 // reduce system to a polynomial from which all solution are extracted 931 // by doing that a lot of global background variables are set 932 TY_ReduceNonlinearSystem( Z1, Z2, K1, K2, phi, debug ); 909 933 910 934 // vector of real and imaginary coefficients in order of decreasing powers 911 int i;912 935 for ( i = 0; i <= degree; i++ ) 913 936 { … … 921 944 922 945 // show the result if in debug mode 923 double x, y; 946 /* 924 947 if ( debug ) 925 948 { … … 937 960 XOPNotice(buf); 938 961 } 962 */ 939 963 940 964 // select real roots and those satisfying Q(x) != 0 and W(x) != 0 … … 947 971 // and g(r) of the correct root should have the minimum average value 948 972 // inside the hardcore 949 double var_a, var_b, var_c1, var_c2, var_d1, var_d2; 950 double sol_a[22], sol_b[22], sol_c1[22], sol_c2[22], sol_d1[22], sol_d2[22]; 951 952 int j = 0; 973 953 974 for ( i = 0; i < degree; i++ ) 954 975 { … … 976 997 sol_d1[j] = var_d1; 977 998 sol_d2[j] = var_d2; 978 999 /* 979 1000 if ( debug ) 980 1001 { … … 991 1012 XOPNotice(buf); 992 1013 } 1014 */ 993 1015 j++; 994 1016 } … … 996 1018 } 997 1019 // number remaining roots 998 int n_roots = j; 1020 n_roots = j; 1021 999 1022 1000 1023 // if there is still more than one root left, than choose the one with the minimum … … 1004 1027 // the number of q values should be a power of 2 1005 1028 // in order to speed up the FFT 1006 intn = 1 << 14;1029 n = 1 << 14; 1007 1030 1008 1031 // the maximum q value should be large enough 1009 1032 // to enable a reasoble approximation of g(r) 1010 doubleqmax = 16 * 10 * 2 * M_PI;1011 d ouble q, dq = qmax / ( n - 1 );1033 qmax = 16 * 10 * 2 * M_PI; 1034 dq = qmax / ( n - 1 ); 1012 1035 1013 // step size for g(r) 1014 double dr; 1036 // step size for g(r) = dr 1015 1037 1016 1038 // allocate memory for pair correlation function g(r) 1017 1039 // and structure factor S(q) 1018 double*sq = malloc( sizeof( double ) * n );1019 double*gr = malloc( sizeof( double ) * n );1040 sq = malloc( sizeof( double ) * n ); 1041 gr = malloc( sizeof( double ) * n ); 1020 1042 1021 1043 // loop over all remaining roots 1022 doublemin = INFINITY;1023 intselected_root = 10;1024 doublesum = 0;1044 min = INFINITY; 1045 selected_root = 10; 1046 sum = 0; 1025 1047 for ( j = 0; j < n_roots; j++) 1026 1048 { … … 1030 1052 q = dq * i; 1031 1053 sq[i] = SqTwoYukawa( q, Z1, Z2, K1, K2, phi, sol_a[j], sol_b[j], sol_c1[j], sol_c2[j], sol_d1[j], sol_d2[j] ); 1032 1054 /* 1033 1055 if(i<20 && debug) { 1034 1056 sprintf(buf, "after SqTwoYukawa: s(q) = %g\r",sq[i] ); 1035 1057 XOPNotice(buf); 1036 1058 } 1059 */ 1037 1060 } 1038 1061 … … 1046 1069 for (i = 0; i < floor( 1. / dr ); i++ ) { 1047 1070 sum += fabs( gr[i] ); 1048 1071 /* 1049 1072 if(i<20 && debug) { 1050 1073 sprintf(buf, "g(r) in core = %g\r",fabs(gr[i])); 1051 1074 XOPNotice(buf); 1052 1075 } 1053 1076 */ 1054 1077 } 1055 1078
Note: See TracChangeset
for help on using the changeset viewer.