Ignore:
Timestamp:
Nov 18, 2008 3:26:32 PM (14 years ago)
Author:
srkline
Message:

Additions to the library of the 2008 model functions. Direct proting of the Igor code, duplicated by these XOPs

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sans/XOP_Dev/SANSAnalysis/lib/libCylinder.c

    r356 r453  
    243243         
    244244        Pi = 4.0*atan(1.0); 
    245         va = 0; 
    246         vb = 1;         //orintational average, outer integral 
    247         vaj = 0; 
    248         vbj = 1;                //endpoints of inner integral 
     245        va = 0.0; 
     246        vb = 1.0;               //orintational average, outer integral 
     247        vaj = 0.0; 
     248        vbj = 1.0;              //endpoints of inner integral 
    249249         
    250250        summ = 0.0;                     //initialize intergral 
     
    260260        for(i=0;i<nordi;i++) { 
    261261                //setup inner integral over the ellipsoidal cross-section 
    262                 summj=0; 
     262                summj=0.0; 
    263263                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy 
    264264                for(j=0;j<nordj;j++) { 
     
    280280        answer *= delrho*delrho; 
    281281        //normalize by ellipsoid volume 
    282         answer *= 4*Pi/3*aa*bb*cc; 
     282        answer *= 4.0*Pi/3.0*aa*bb*cc; 
    283283        //convert to [cm-1] 
    284284        answer *= 1.0e8; 
     
    22192219        return(ans); 
    22202220} 
     2221 
     2222/*      Lamellar_ParaCrystal - Pedersen's model 
     2223 
     2224*/ 
     2225double 
     2226Lamellar_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.0-xn)*Snq; 
     2274         
     2275        //and the independent contribution 
     2276        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg)); 
     2277         
     2278        //the limit when NL approaches infinity 
     2279//      Zq = (1-ww^2)/(1+ww^2-2*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 
     2291double 
     2292paraCryst_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*ww-2.0*ww*cos(qval*davg)),2) ); 
     2297         
     2298        return(Snq); 
     2299} 
     2300 
     2301 
     2302double 
     2303paraCryst_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)(nl-1)*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 
     2318Uses 76 pt Gaussian quadrature for both integrals 
     2319*/ 
     2320double 
     2321Spherocylinder(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 = sldc-slds; 
     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 cross-section 
     2365                summj=0.0; 
     2366                zi = ( Gauss76Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/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// 
     2408double 
     2409SphCyl_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.0-tt*tt); 
     2426         
     2427        val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; 
     2428         
     2429        return(val); 
     2430} 
     2431 
     2432 
     2433/*      Convex Lens  : special case where L ~ 0 and hDist < 0 
     2434 
     2435Uses 76 pt Gaussian quadrature for both integrals 
     2436*/ 
     2437double 
     2438ConvexLens(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*endRad-rad*rad));         //by definition for this model 
     2471         
     2472        contr = sldc-slds; 
     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 cross-section 
     2484                summj=0.0; 
     2485                zi = ( Gauss76Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/2.0*summ; 
     2514 
     2515        hh = fabs(hDist);               //need positive value for spherical cap volume 
     2516        answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(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  
     2529Uses 76 pt Gaussian quadrature for both integrals 
     2530*/ 
     2531double 
     2532CappedCylinder(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*endRad-rad*rad));         //by definition for this model 
     2555         
     2556        contr = sldc-slds; 
     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 cross-section 
     2568                summj=0.0; 
     2569                zi = ( Gauss76Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/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*(endRad-hh)*(endRad-hh)*(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// 
     2613double 
     2614ConvLens_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*endRad-rad*rad)); 
     2628                 
     2629        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 
     2630        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 
     2631         
     2632        val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; 
     2633         
     2634        return(val); 
     2635} 
     2636 
     2637 
     2638/*      Dumbbell  : special case where L ~ 0 and hDist > 0 
     2639 
     2640Uses 76 pt Gaussian quadrature for both integrals 
     2641*/ 
     2642double 
     2643Dumbbell(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*endRad-rad*rad));              //by definition for this model 
     2676         
     2677        contr = sldc-slds; 
     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 cross-section 
     2689                summj=0.0; 
     2690                zi = ( Gauss76Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/2.0*summ; 
     2719 
     2720        answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*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 
     2734Uses 76 pt Gaussian quadrature for both integrals 
     2735*/ 
     2736double 
     2737Barbell(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*endRad-rad*rad));              //by definition for this model 
     2760         
     2761        contr = sldc-slds; 
     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 cross-section 
     2773                summj=0.0; 
     2774                zi = ( Gauss76Z[i]*(vb-va) + 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]*(vbj-vaj) + 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 = (vbj-vaj)/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 = (vb-va)/2.0*summ; 
     2803 
     2804        answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*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// 
     2819double 
     2820Dumb_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*endRad-rad*rad)); 
     2834                 
     2835        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0); 
     2836        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 
     2837         
     2838        val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; 
     2839         
     2840        return(val); 
     2841} 
Note: See TracChangeset for help on using the changeset viewer.