Changeset 453 for sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
 Timestamp:
 Nov 18, 2008 3:26:32 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c
r356 r453 243 243 244 244 Pi = 4.0*atan(1.0); 245 va = 0 ;246 vb = 1 ; //orintational average, outer integral247 vaj = 0 ;248 vbj = 1 ; //endpoints of inner integral245 va = 0.0; 246 vb = 1.0; //orintational average, outer integral 247 vaj = 0.0; 248 vbj = 1.0; //endpoints of inner integral 249 249 250 250 summ = 0.0; //initialize intergral … … 260 260 for(i=0;i<nordi;i++) { 261 261 //setup inner integral over the ellipsoidal crosssection 262 summj=0 ;262 summj=0.0; 263 263 zi = ( Gauss76Z[i]*(vbva) + va + vb )/2.0; //the "x" dummy 264 264 for(j=0;j<nordj;j++) { … … 280 280 answer *= delrho*delrho; 281 281 //normalize by ellipsoid volume 282 answer *= 4 *Pi/3*aa*bb*cc;282 answer *= 4.0*Pi/3.0*aa*bb*cc; 283 283 //convert to [cm1] 284 284 answer *= 1.0e8; … … 2219 2219 return(ans); 2220 2220 } 2221 2222 /* Lamellar_ParaCrystal  Pedersen's model 2223 2224 */ 2225 double 2226 Lamellar_ParaCrystal(double w[], double q) 2227 { 2228 // Input (fitting) variables are: 2229 //[0] scale factor 2230 //[1] thickness 2231 //[2] number of layers 2232 //[3] spacing between layers 2233 //[4] polydispersity of spacing 2234 //[5] SLD lamellar 2235 //[6] SLD solvent 2236 //[7] incoherent background 2237 // give them nice names 2238 double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn; 2239 double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi; 2240 long n1,n2; 2241 2242 pi = 4.0*atan(1.0); 2243 scale = w[0]; 2244 th = w[1]; 2245 nl = w[2]; 2246 davg = w[3]; 2247 pd = w[4]; 2248 sldLayer = w[5]; 2249 sldSolvent = w[6]; 2250 bkg = w[7]; 2251 2252 contr = w[5]  w[6]; 2253 qval = q; 2254 2255 //get the fractional part of nl, to determine the "mixing" of N's 2256 2257 n1 = trunc(nl); //rounds towards zero 2258 n2 = n1 + 1; 2259 xn = (double)n2  nl; //fractional contribution of n1 2260 2261 ww = exp(qval*qval*pd*pd*davg*davg/2.0); 2262 2263 //calculate the n1 contribution 2264 an = paraCryst_an(ww,qval,davg,n1); 2265 Snq = paraCryst_sn(ww,qval,davg,n1,an); 2266 2267 Znq = xn*Snq; 2268 2269 //calculate the n2 contribution 2270 an = paraCryst_an(ww,qval,davg,n2); 2271 Snq = paraCryst_sn(ww,qval,davg,n2,an); 2272 2273 Znq += (1.0xn)*Snq; 2274 2275 //and the independent contribution 2276 Znq += (1.0ww*ww)/(1.0+ww*ww2.0*ww*cos(qval*davg)); 2277 2278 //the limit when NL approaches infinity 2279 // Zq = (1ww^2)/(1+ww^22*ww*cos(qval*davg)) 2280 2281 xi = th/2.0; //use 1/2 the bilayer thickness 2282 Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi)); 2283 2284 inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval); 2285 inten *= 1.0e8; 2286 2287 return(scale*inten+bkg); 2288 } 2289 2290 // functions for the lamellar paracrystal model 2291 double 2292 paraCryst_sn(double ww, double qval, double davg, long nl, double an) { 2293 2294 double Snq; 2295 2296 Snq = an/( (double)nl*pow((1.0+ww*ww2.0*ww*cos(qval*davg)),2) ); 2297 2298 return(Snq); 2299 } 2300 2301 2302 double 2303 paraCryst_an(double ww, double qval, double davg, long nl) { 2304 2305 double an; 2306 2307 an = 4.0*ww*ww  2.0*(ww*ww*ww+ww)*cos(qval*davg); 2308 an = 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg); 2309 an += 2.0*pow(ww,(nl+3))*cos((double)(nl1)*qval*davg); 2310 an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg); 2311 2312 return(an); 2313 } 2314 2315 2316 /* Spherocylinder : 2317 2318 Uses 76 pt Gaussian quadrature for both integrals 2319 */ 2320 double 2321 Spherocylinder(double w[], double x) 2322 { 2323 int i,j; 2324 double Pi; 2325 double scale,contr,bkg,sldc,slds; 2326 double len,rad,hDist,endRad; 2327 int nordi=76; //order of integration 2328 int nordj=76; 2329 double va,vb; //upper and lower integration limits 2330 double summ,zi,yyy,answer; //running tally of integration 2331 double summj,vaj,vbj,zij; //for the inner integration 2332 double SphCyl_tmp[7],arg1,arg2,inner; 2333 2334 2335 scale = w[0]; 2336 rad = w[1]; 2337 len = w[2]; 2338 sldc = w[3]; 2339 slds = w[4]; 2340 bkg = w[5]; 2341 2342 SphCyl_tmp[0] = w[0]; 2343 SphCyl_tmp[1] = w[1]; 2344 SphCyl_tmp[2] = w[2]; 2345 SphCyl_tmp[3] = w[1]; //end radius is same as cylinder radius 2346 SphCyl_tmp[4] = w[3]; 2347 SphCyl_tmp[5] = w[4]; 2348 SphCyl_tmp[6] = w[5]; 2349 2350 hDist = 0; //by definition for this model 2351 endRad = rad; 2352 2353 contr = sldcslds; 2354 2355 Pi = 4.0*atan(1.0); 2356 va = 0.0; 2357 vb = Pi/2.0; //orintational average, outer integral 2358 vaj = 1.0*hDist/endRad; 2359 vbj = 1.0; //endpoints of inner integral 2360 2361 summ = 0.0; //initialize intergral 2362 2363 for(i=0;i<nordi;i++) { 2364 //setup inner integral over the ellipsoidal crosssection 2365 summj=0.0; 2366 zi = ( Gauss76Z[i]*(vbva) + va + vb )/2.0; //the "theta" dummy 2367 2368 for(j=0;j<nordj;j++) { 2369 //20 gauss points for the inner integral 2370 zij = ( Gauss76Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the "t" dummy 2371 yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi); 2372 summj += yyy; 2373 } 2374 //now calculate the value of the inner integral 2375 inner = (vbjvaj)/2.0*summj; 2376 inner *= 4.0*Pi*endRad*endRad*endRad; 2377 2378 //now calculate outer integrand 2379 arg1 = x*len/2.0*cos(zi); 2380 arg2 = x*rad*sin(zi); 2381 yyy = inner; 2382 2383 if(arg1 == 0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2384 yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 2385 } else { 2386 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 2387 } 2388 yyy *= yyy; 2389 yyy *= sin(zi); // = A(q)^2*sin(theta) 2390 yyy *= Gauss76Wt[i]; 2391 summ += yyy; 2392 } //final scaling is done at the end of the function, after the NT_FP64 case 2393 2394 answer = (vbva)/2.0*summ; 2395 2396 answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0; //divide by volume 2397 answer *= 1.0e8; //convert to cm^1 2398 answer *= contr*contr; 2399 answer *= scale; 2400 answer += bkg; 2401 2402 return answer; 2403 } 2404 2405 2406 // inner integral of the sphereocylinder model, special case of hDist = 0 2407 // 2408 double 2409 SphCyl_kernel(double w[], double x, double tt, double theta) { 2410 2411 double val,arg1,arg2; 2412 double scale,bkg,sldc,slds; 2413 double len,rad,hDist,endRad; 2414 scale = w[0]; 2415 rad = w[1]; 2416 len = w[2]; 2417 endRad = w[3]; 2418 sldc = w[4]; 2419 slds = w[5]; 2420 bkg = w[6]; 2421 2422 hDist = 0.0; 2423 2424 arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 2425 arg2 = x*endRad*sin(theta)*sqrt(1.0tt*tt); 2426 2427 val = cos(arg1)*(1.0tt*tt)*NR_BessJ1(arg2)/arg2; 2428 2429 return(val); 2430 } 2431 2432 2433 /* Convex Lens : special case where L ~ 0 and hDist < 0 2434 2435 Uses 76 pt Gaussian quadrature for both integrals 2436 */ 2437 double 2438 ConvexLens(double w[], double x) 2439 { 2440 int i,j; 2441 double Pi; 2442 double scale,contr,bkg,sldc,slds; 2443 double len,rad,hDist,endRad; 2444 int nordi=76; //order of integration 2445 int nordj=76; 2446 double va,vb; //upper and lower integration limits 2447 double summ,zi,yyy,answer; //running tally of integration 2448 double summj,vaj,vbj,zij; //for the inner integration 2449 double CLens_tmp[7],arg1,arg2,inner,hh; 2450 2451 2452 scale = w[0]; 2453 rad = w[1]; 2454 // len = w[2] 2455 endRad = w[2]; 2456 sldc = w[3]; 2457 slds = w[4]; 2458 bkg = w[5]; 2459 2460 len = 0.01; 2461 2462 CLens_tmp[0] = w[0]; 2463 CLens_tmp[1] = w[1]; 2464 CLens_tmp[2] = len; //length is some small number, essentially zero 2465 CLens_tmp[3] = w[2]; 2466 CLens_tmp[4] = w[3]; 2467 CLens_tmp[5] = w[4]; 2468 CLens_tmp[6] = w[5]; 2469 2470 hDist = 1.0*sqrt(fabs(endRad*endRadrad*rad)); //by definition for this model 2471 2472 contr = sldcslds; 2473 2474 Pi = 4.0*atan(1.0); 2475 va = 0.0; 2476 vb = Pi/2.0; //orintational average, outer integral 2477 vaj = 1.0*hDist/endRad; 2478 vbj = 1.0; //endpoints of inner integral 2479 2480 summ = 0.0; //initialize intergral 2481 2482 for(i=0;i<nordi;i++) { 2483 //setup inner integral over the ellipsoidal crosssection 2484 summj=0.0; 2485 zi = ( Gauss76Z[i]*(vbva) + va + vb )/2.0; //the "theta" dummy 2486 2487 for(j=0;j<nordj;j++) { 2488 //20 gauss points for the inner integral 2489 zij = ( Gauss76Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the "t" dummy 2490 yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi); 2491 summj += yyy; 2492 } 2493 //now calculate the value of the inner integral 2494 inner = (vbjvaj)/2.0*summj; 2495 inner *= 4.0*Pi*endRad*endRad*endRad; 2496 2497 //now calculate outer integrand 2498 arg1 = x*len/2.0*cos(zi); 2499 arg2 = x*rad*sin(zi); 2500 yyy = inner; 2501 2502 if(arg1 == 0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2503 yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 2504 } else { 2505 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 2506 } 2507 yyy *= yyy; 2508 yyy *= sin(zi); // = A(q)^2*sin(theta) 2509 yyy *= Gauss76Wt[i]; 2510 summ += yyy; 2511 } //final scaling is done at the end of the function, after the NT_FP64 case 2512 2513 answer = (vbva)/2.0*summ; 2514 2515 hh = fabs(hDist); //need positive value for spherical cap volume 2516 answer /= 2.0*(1.0/3.0*Pi*(endRadhh)*(endRadhh)*(2.0*endRad+hh)); //divide by volume 2517 answer *= 1.0e8; //convert to cm^1 2518 answer *= contr*contr; 2519 answer *= scale; 2520 answer += bkg; 2521 2522 return answer; 2523 } 2524 2525 /* Capped Cylinder : special case where L is nonzero and hDist < 0 2526 2527  uses the same Kernel as the Convex Lens 2528 2529 Uses 76 pt Gaussian quadrature for both integrals 2530 */ 2531 double 2532 CappedCylinder(double w[], double x) 2533 { 2534 int i,j; 2535 double Pi; 2536 double scale,contr,bkg,sldc,slds; 2537 double len,rad,hDist,endRad; 2538 int nordi=76; //order of integration 2539 int nordj=76; 2540 double va,vb; //upper and lower integration limits 2541 double summ,zi,yyy,answer; //running tally of integration 2542 double summj,vaj,vbj,zij; //for the inner integration 2543 double arg1,arg2,inner,hh; 2544 2545 2546 scale = w[0]; 2547 rad = w[1]; 2548 len = w[2]; 2549 endRad = w[3]; 2550 sldc = w[4]; 2551 slds = w[5]; 2552 bkg = w[6]; 2553 2554 hDist = 1.0*sqrt(fabs(endRad*endRadrad*rad)); //by definition for this model 2555 2556 contr = sldcslds; 2557 2558 Pi = 4.0*atan(1.0); 2559 va = 0.0; 2560 vb = Pi/2.0; //orintational average, outer integral 2561 vaj = 1.0*hDist/endRad; 2562 vbj = 1.0; //endpoints of inner integral 2563 2564 summ = 0.0; //initialize intergral 2565 2566 for(i=0;i<nordi;i++) { 2567 //setup inner integral over the ellipsoidal crosssection 2568 summj=0.0; 2569 zi = ( Gauss76Z[i]*(vbva) + va + vb )/2.0; //the "theta" dummy 2570 2571 for(j=0;j<nordj;j++) { 2572 //20 gauss points for the inner integral 2573 zij = ( Gauss76Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the "t" dummy 2574 yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi); //uses the same kernel as ConvexLens, except here L != 0 2575 summj += yyy; 2576 } 2577 //now calculate the value of the inner integral 2578 inner = (vbjvaj)/2.0*summj; 2579 inner *= 4.0*Pi*endRad*endRad*endRad; 2580 2581 //now calculate outer integrand 2582 arg1 = x*len/2.0*cos(zi); 2583 arg2 = x*rad*sin(zi); 2584 yyy = inner; 2585 2586 if(arg1 == 0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2587 yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 2588 } else { 2589 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 2590 } 2591 yyy *= yyy; 2592 yyy *= sin(zi); // = A(q)^2*sin(theta) 2593 yyy *= Gauss76Wt[i]; 2594 summ += yyy; 2595 } //final scaling is done at the end of the function, after the NT_FP64 case 2596 2597 answer = (vbva)/2.0*summ; 2598 2599 hh = fabs(hDist); //need positive value for spherical cap volume 2600 answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRadhh)*(endRadhh)*(2.0*endRad+hh)); //divide by volume 2601 answer *= 1.0e8; //convert to cm^1 2602 answer *= contr*contr; 2603 answer *= scale; 2604 answer += bkg; 2605 2606 return answer; 2607 } 2608 2609 2610 2611 // inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0 2612 // 2613 double 2614 ConvLens_kernel(double w[], double x, double tt, double theta) { 2615 2616 double val,arg1,arg2; 2617 double scale,bkg,sldc,slds; 2618 double len,rad,hDist,endRad; 2619 scale = w[0]; 2620 rad = w[1]; 2621 len = w[2]; 2622 endRad = w[3]; 2623 sldc = w[4]; 2624 slds = w[5]; 2625 bkg = w[6]; 2626 2627 hDist = 1.0*sqrt(fabs(endRad*endRadrad*rad)); 2628 2629 arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 2630 arg2 = x*endRad*sin(theta)*sqrt(1.0tt*tt); 2631 2632 val = cos(arg1)*(1.0tt*tt)*NR_BessJ1(arg2)/arg2; 2633 2634 return(val); 2635 } 2636 2637 2638 /* Dumbbell : special case where L ~ 0 and hDist > 0 2639 2640 Uses 76 pt Gaussian quadrature for both integrals 2641 */ 2642 double 2643 Dumbbell(double w[], double x) 2644 { 2645 int i,j; 2646 double Pi; 2647 double scale,contr,bkg,sldc,slds; 2648 double len,rad,hDist,endRad; 2649 int nordi=76; //order of integration 2650 int nordj=76; 2651 double va,vb; //upper and lower integration limits 2652 double summ,zi,yyy,answer; //running tally of integration 2653 double summj,vaj,vbj,zij; //for the inner integration 2654 double Dumb_tmp[7],arg1,arg2,inner; 2655 2656 2657 scale = w[0]; 2658 rad = w[1]; 2659 // len = w[2] 2660 endRad = w[2]; 2661 sldc = w[3]; 2662 slds = w[4]; 2663 bkg = w[5]; 2664 2665 len = 0.01; 2666 2667 Dumb_tmp[0] = w[0]; 2668 Dumb_tmp[1] = w[1]; 2669 Dumb_tmp[2] = len; //length is some small number, essentially zero 2670 Dumb_tmp[3] = w[2]; 2671 Dumb_tmp[4] = w[3]; 2672 Dumb_tmp[5] = w[4]; 2673 Dumb_tmp[6] = w[5]; 2674 2675 hDist = sqrt(fabs(endRad*endRadrad*rad)); //by definition for this model 2676 2677 contr = sldcslds; 2678 2679 Pi = 4.0*atan(1.0); 2680 va = 0.0; 2681 vb = Pi/2.0; //orintational average, outer integral 2682 vaj = 1.0*hDist/endRad; 2683 vbj = 1.0; //endpoints of inner integral 2684 2685 summ = 0.0; //initialize intergral 2686 2687 for(i=0;i<nordi;i++) { 2688 //setup inner integral over the ellipsoidal crosssection 2689 summj=0.0; 2690 zi = ( Gauss76Z[i]*(vbva) + va + vb )/2.0; //the "theta" dummy 2691 2692 for(j=0;j<nordj;j++) { 2693 //20 gauss points for the inner integral 2694 zij = ( Gauss76Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the "t" dummy 2695 yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi); 2696 summj += yyy; 2697 } 2698 //now calculate the value of the inner integral 2699 inner = (vbjvaj)/2.0*summj; 2700 inner *= 4.0*Pi*endRad*endRad*endRad; 2701 2702 //now calculate outer integrand 2703 arg1 = x*len/2.0*cos(zi); 2704 arg2 = x*rad*sin(zi); 2705 yyy = inner; 2706 2707 if(arg1 == 0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2708 yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 2709 } else { 2710 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 2711 } 2712 yyy *= yyy; 2713 yyy *= sin(zi); // = A(q)^2*sin(theta) 2714 yyy *= Gauss76Wt[i]; 2715 summ += yyy; 2716 } //final scaling is done at the end of the function, after the NT_FP64 case 2717 2718 answer = (vbva)/2.0*summ; 2719 2720 answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDisthDist*hDist*hDist/3.0); //divide by volume 2721 answer *= 1.0e8; //convert to cm^1 2722 answer *= contr*contr; 2723 answer *= scale; 2724 answer += bkg; 2725 2726 return answer; 2727 } 2728 2729 2730 /* Barbell : "normal" case where L is nonzero 0 and hDist > 0 2731 2732  uses the same kernel as the Dumbbell case 2733 2734 Uses 76 pt Gaussian quadrature for both integrals 2735 */ 2736 double 2737 Barbell(double w[], double x) 2738 { 2739 int i,j; 2740 double Pi; 2741 double scale,contr,bkg,sldc,slds; 2742 double len,rad,hDist,endRad; 2743 int nordi=76; //order of integration 2744 int nordj=76; 2745 double va,vb; //upper and lower integration limits 2746 double summ,zi,yyy,answer; //running tally of integration 2747 double summj,vaj,vbj,zij; //for the inner integration 2748 double arg1,arg2,inner; 2749 2750 2751 scale = w[0]; 2752 rad = w[1]; 2753 len = w[2]; 2754 endRad = w[3]; 2755 sldc = w[4]; 2756 slds = w[5]; 2757 bkg = w[6]; 2758 2759 hDist = sqrt(fabs(endRad*endRadrad*rad)); //by definition for this model 2760 2761 contr = sldcslds; 2762 2763 Pi = 4.0*atan(1.0); 2764 va = 0.0; 2765 vb = Pi/2.0; //orintational average, outer integral 2766 vaj = 1.0*hDist/endRad; 2767 vbj = 1.0; //endpoints of inner integral 2768 2769 summ = 0.0; //initialize intergral 2770 2771 for(i=0;i<nordi;i++) { 2772 //setup inner integral over the ellipsoidal crosssection 2773 summj=0.0; 2774 zi = ( Gauss76Z[i]*(vbva) + va + vb )/2.0; //the "theta" dummy 2775 2776 for(j=0;j<nordj;j++) { 2777 //20 gauss points for the inner integral 2778 zij = ( Gauss76Z[j]*(vbjvaj) + vaj + vbj )/2.0; //the "t" dummy 2779 yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi); //uses the same Kernel as the Dumbbell, here L>0 2780 summj += yyy; 2781 } 2782 //now calculate the value of the inner integral 2783 inner = (vbjvaj)/2.0*summj; 2784 inner *= 4.0*Pi*endRad*endRad*endRad; 2785 2786 //now calculate outer integrand 2787 arg1 = x*len/2.0*cos(zi); 2788 arg2 = x*rad*sin(zi); 2789 yyy = inner; 2790 2791 if(arg1 == 0) { //limiting value of sinc(0) is 1; sinc is not defined in math.h 2792 yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 2793 } else { 2794 yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 2795 } 2796 yyy *= yyy; 2797 yyy *= sin(zi); // = A(q)^2*sin(theta) 2798 yyy *= Gauss76Wt[i]; 2799 summ += yyy; 2800 } //final scaling is done at the end of the function, after the NT_FP64 case 2801 2802 answer = (vbva)/2.0*summ; 2803 2804 answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDisthDist*hDist*hDist/3.0); //divide by volume 2805 answer *= 1.0e8; //convert to cm^1 2806 answer *= contr*contr; 2807 answer *= scale; 2808 answer += bkg; 2809 2810 return answer; 2811 } 2812 2813 2814 2815 // inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0 2816 // 2817 // inner integral of the Barbell model if L is nonzero 2818 // 2819 double 2820 Dumb_kernel(double w[], double x, double tt, double theta) { 2821 2822 double val,arg1,arg2; 2823 double scale,bkg,sldc,slds; 2824 double len,rad,hDist,endRad; 2825 scale = w[0]; 2826 rad = w[1]; 2827 len = w[2]; 2828 endRad = w[3]; 2829 sldc = w[4]; 2830 slds = w[5]; 2831 bkg = w[6]; 2832 2833 hDist = sqrt(fabs(endRad*endRadrad*rad)); 2834 2835 arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 2836 arg2 = x*endRad*sin(theta)*sqrt(1.0tt*tt); 2837 2838 val = cos(arg1)*(1.0tt*tt)*NR_BessJ1(arg2)/arg2; 2839 2840 return(val); 2841 }
Note: See TracChangeset
for help on using the changeset viewer.