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

Location:
sans/XOP_Dev/SANSAnalysis/lib
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • sans/XOP_Dev/SANSAnalysis/lib/GaussWeights.h

    r97 r453  
    212212}; 
    213213 
     214static double Gauss150Z[150]={ 
     215        -0.9998723404457334, 
     216        -0.9993274305065947, 
     217        -0.9983473449340834, 
     218        -0.9969322929775997, 
     219        -0.9950828645255290, 
     220        -0.9927998590434373, 
     221        -0.9900842691660192, 
     222        -0.9869372772712794, 
     223        -0.9833602541697529, 
     224        -0.9793547582425894, 
     225        -0.9749225346595943, 
     226        -0.9700655145738374, 
     227        -0.9647858142586956, 
     228        -0.9590857341746905, 
     229        -0.9529677579610971, 
     230        -0.9464345513503147, 
     231        -0.9394889610042837, 
     232        -0.9321340132728527, 
     233        -0.9243729128743136, 
     234        -0.9162090414984952, 
     235        -0.9076459563329236, 
     236        -0.8986873885126239, 
     237        -0.8893372414942055, 
     238        -0.8795995893549102, 
     239        -0.8694786750173527, 
     240        -0.8589789084007133, 
     241        -0.8481048644991847, 
     242        -0.8368612813885015, 
     243        -0.8252530581614230, 
     244        -0.8132852527930605, 
     245        -0.8009630799369827, 
     246        -0.7882919086530552, 
     247        -0.7752772600680049, 
     248        -0.7619248049697269, 
     249        -0.7482403613363824, 
     250        -0.7342298918013638, 
     251        -0.7198995010552305, 
     252        -0.7052554331857488, 
     253        -0.6903040689571928, 
     254        -0.6750519230300931, 
     255        -0.6595056411226444, 
     256        -0.6436719971150083, 
     257        -0.6275578900977726, 
     258        -0.6111703413658551, 
     259        -0.5945164913591590, 
     260        -0.5776035965513142, 
     261        -0.5604390262878617, 
     262        -0.5430302595752546, 
     263        -0.5253848818220803, 
     264        -0.5075105815339176, 
     265        -0.4894151469632753, 
     266        -0.4711064627160663, 
     267        -0.4525925063160997, 
     268        -0.4338813447290861, 
     269        -0.4149811308476706, 
     270        -0.3959000999390257, 
     271        -0.3766465660565522, 
     272        -0.3572289184172501, 
     273        -0.3376556177463400, 
     274        -0.3179351925907259, 
     275        -0.2980762356029071, 
     276        -0.2780873997969574, 
     277        -0.2579773947782034, 
     278        -0.2377549829482451, 
     279        -0.2174289756869712, 
     280        -0.1970082295132342, 
     281        -0.1765016422258567, 
     282        -0.1559181490266516, 
     283        -0.1352667186271445, 
     284        -0.1145563493406956, 
     285        -0.0937960651617229, 
     286        -0.0729949118337358, 
     287        -0.0521619529078925, 
     288        -0.0313062657937972, 
     289        -0.0104369378042598, 
     290        0.0104369378042598, 
     291        0.0313062657937972, 
     292        0.0521619529078925, 
     293        0.0729949118337358, 
     294        0.0937960651617229, 
     295        0.1145563493406956, 
     296        0.1352667186271445, 
     297        0.1559181490266516, 
     298        0.1765016422258567, 
     299        0.1970082295132342, 
     300        0.2174289756869712, 
     301        0.2377549829482451, 
     302        0.2579773947782034, 
     303        0.2780873997969574, 
     304        0.2980762356029071, 
     305        0.3179351925907259, 
     306        0.3376556177463400, 
     307        0.3572289184172501, 
     308        0.3766465660565522, 
     309        0.3959000999390257, 
     310        0.4149811308476706, 
     311        0.4338813447290861, 
     312        0.4525925063160997, 
     313        0.4711064627160663, 
     314        0.4894151469632753, 
     315        0.5075105815339176, 
     316        0.5253848818220803, 
     317        0.5430302595752546, 
     318        0.5604390262878617, 
     319        0.5776035965513142, 
     320        0.5945164913591590, 
     321        0.6111703413658551, 
     322        0.6275578900977726, 
     323        0.6436719971150083, 
     324        0.6595056411226444, 
     325        0.6750519230300931, 
     326        0.6903040689571928, 
     327        0.7052554331857488, 
     328        0.7198995010552305, 
     329        0.7342298918013638, 
     330        0.7482403613363824, 
     331        0.7619248049697269, 
     332        0.7752772600680049, 
     333        0.7882919086530552, 
     334        0.8009630799369827, 
     335        0.8132852527930605, 
     336        0.8252530581614230, 
     337        0.8368612813885015, 
     338        0.8481048644991847, 
     339        0.8589789084007133, 
     340        0.8694786750173527, 
     341        0.8795995893549102, 
     342        0.8893372414942055, 
     343        0.8986873885126239, 
     344        0.9076459563329236, 
     345        0.9162090414984952, 
     346        0.9243729128743136, 
     347        0.9321340132728527, 
     348        0.9394889610042837, 
     349        0.9464345513503147, 
     350        0.9529677579610971, 
     351        0.9590857341746905, 
     352        0.9647858142586956, 
     353        0.9700655145738374, 
     354        0.9749225346595943, 
     355        0.9793547582425894, 
     356        0.9833602541697529, 
     357        0.9869372772712794, 
     358        0.9900842691660192, 
     359        0.9927998590434373, 
     360        0.9950828645255290, 
     361        0.9969322929775997, 
     362        0.9983473449340834, 
     363        0.9993274305065947, 
     364        0.9998723404457334 
     365}; 
    214366 
     367static double Gauss150Wt[150]={ 
     368        0.0003276086705538, 
     369        0.0007624720924706, 
     370        0.0011976474864367, 
     371        0.0016323569986067, 
     372        0.0020663664924131, 
     373        0.0024994789888943, 
     374        0.0029315036836558, 
     375        0.0033622516236779, 
     376        0.0037915348363451, 
     377        0.0042191661429919, 
     378        0.0046449591497966, 
     379        0.0050687282939456, 
     380        0.0054902889094487, 
     381        0.0059094573005900, 
     382        0.0063260508184704, 
     383        0.0067398879387430, 
     384        0.0071507883396855, 
     385        0.0075585729801782, 
     386        0.0079630641773633, 
     387        0.0083640856838475, 
     388        0.0087614627643580, 
     389        0.0091550222717888, 
     390        0.0095445927225849, 
     391        0.0099300043714212, 
     392        0.0103110892851360, 
     393        0.0106876814158841, 
     394        0.0110596166734735, 
     395        0.0114267329968529, 
     396        0.0117888704247183, 
     397        0.0121458711652067, 
     398        0.0124975796646449, 
     399        0.0128438426753249, 
     400        0.0131845093222756, 
     401        0.0135194311690004, 
     402        0.0138484622795371, 
     403        0.0141714592928592, 
     404        0.0144882814685445, 
     405        0.0147987907597169, 
     406        0.0151028518701744, 
     407        0.0154003323133401, 
     408        0.0156911024699895, 
     409        0.0159750356447283, 
     410        0.0162520081211971, 
     411        0.0165218992159766, 
     412        0.0167845913311726, 
     413        0.0170399700056559, 
     414        0.0172879239649355, 
     415        0.0175283451696437, 
     416        0.0177611288626114, 
     417        0.0179861736145128, 
     418        0.0182033813680609, 
     419        0.0184126574807331, 
     420        0.0186139107660094, 
     421        0.0188070535331042, 
     422        0.0189920016251754, 
     423        0.0191686744559934, 
     424        0.0193369950450545, 
     425        0.0194968900511231, 
     426        0.0196482898041878, 
     427        0.0197911283358190, 
     428        0.0199253434079123, 
     429        0.0200508765398072, 
     430        0.0201676730337687, 
     431        0.0202756819988200, 
     432        0.0203748563729175, 
     433        0.0204651529434560, 
     434        0.0205465323660984, 
     435        0.0206189591819181, 
     436        0.0206824018328499, 
     437        0.0207368326754401, 
     438        0.0207822279928917, 
     439        0.0208185680053983, 
     440        0.0208458368787627, 
     441        0.0208640227312962, 
     442        0.0208731176389954, 
     443        0.0208731176389954, 
     444        0.0208640227312962, 
     445        0.0208458368787627, 
     446        0.0208185680053983, 
     447        0.0207822279928917, 
     448        0.0207368326754401, 
     449        0.0206824018328499, 
     450        0.0206189591819181, 
     451        0.0205465323660984, 
     452        0.0204651529434560, 
     453        0.0203748563729175, 
     454        0.0202756819988200, 
     455        0.0201676730337687, 
     456        0.0200508765398072, 
     457        0.0199253434079123, 
     458        0.0197911283358190, 
     459        0.0196482898041878, 
     460        0.0194968900511231, 
     461        0.0193369950450545, 
     462        0.0191686744559934, 
     463        0.0189920016251754, 
     464        0.0188070535331042, 
     465        0.0186139107660094, 
     466        0.0184126574807331, 
     467        0.0182033813680609, 
     468        0.0179861736145128, 
     469        0.0177611288626114, 
     470        0.0175283451696437, 
     471        0.0172879239649355, 
     472        0.0170399700056559, 
     473        0.0167845913311726, 
     474        0.0165218992159766, 
     475        0.0162520081211971, 
     476        0.0159750356447283, 
     477        0.0156911024699895, 
     478        0.0154003323133401, 
     479        0.0151028518701744, 
     480        0.0147987907597169, 
     481        0.0144882814685445, 
     482        0.0141714592928592, 
     483        0.0138484622795371, 
     484        0.0135194311690004, 
     485        0.0131845093222756, 
     486        0.0128438426753249, 
     487        0.0124975796646449, 
     488        0.0121458711652067, 
     489        0.0117888704247183, 
     490        0.0114267329968529, 
     491        0.0110596166734735, 
     492        0.0106876814158841, 
     493        0.0103110892851360, 
     494        0.0099300043714212, 
     495        0.0095445927225849, 
     496        0.0091550222717888, 
     497        0.0087614627643580, 
     498        0.0083640856838475, 
     499        0.0079630641773633, 
     500        0.0075585729801782, 
     501        0.0071507883396855, 
     502        0.0067398879387430, 
     503        0.0063260508184704, 
     504        0.0059094573005900, 
     505        0.0054902889094487, 
     506        0.0050687282939456, 
     507        0.0046449591497966, 
     508        0.0042191661429919, 
     509        0.0037915348363451, 
     510        0.0033622516236779, 
     511        0.0029315036836558, 
     512        0.0024994789888943, 
     513        0.0020663664924131, 
     514        0.0016323569986067, 
     515        0.0011976474864367, 
     516        0.0007624720924706, 
     517        0.0003276086705538 
     518}; 
    215519 

  • 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} 

  • sans/XOP_Dev/SANSAnalysis/lib/libCylinder.h

    r97 r453  
    2828double LamellarPS(double dp[], double q); 
    2929double LamellarPS_HG(double dp[], double q); 
     30double Lamellar_ParaCrystal(double dp[], double q); 
     31double Spherocylinder(double dp[], double q); 
     32double ConvexLens(double dp[], double q); 
     33double Dumbbell(double dp[], double q); 
     34double CappedCylinder(double dp[], double q); 
     35double Barbell(double dp[], double q); 
     36 
    3037 
    3138/* internal functions */ 
    … …  
    4956double CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length); 
    5057double Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N); 
     58double paraCryst_sn(double ww, double qval, double davg, long nl, double an); 
     59double paraCryst_an(double ww, double qval, double davg, long nl); 
     60double SphCyl_kernel(double w[], double x, double tt, double Theta); 
     61double ConvLens_kernel(double w[], double x, double tt, double theta); 
     62double Dumb_kernel(double w[], double x, double tt, double theta); 
     63 
    5164 
    5265/////////functions for WRC implementation of flexible cylinders 

  • sans/XOP_Dev/SANSAnalysis/lib/libSANSAnalysis.h

    r231 r453  
    3131double LamellarPS(double dp[], double q); 
    3232double LamellarPS_HG(double dp[], double q); 
     33// 
     34double Lamellar_ParaCrystal(double dp[], double q); 
     35double Spherocylinder(double dp[], double q); 
     36double ConvexLens(double dp[], double q); 
     37double Dumbbell(double dp[], double q); 
     38double CappedCylinder(double dp[], double q); 
     39double Barbell(double dp[], double q); 
     40 
    3341 
    3442/* internal functions */ 
    … …  
    7280double BinaryHS_PSF12(double dp[], double q); 
    7381double BinaryHS_PSF22(double dp[], double q); 
     82// 
     83double OneShell(double dp[], double q); 
     84double TwoShell(double dp[], double q); 
     85double ThreeShell(double dp[], double q); 
     86double FourShell(double dp[], double q); 
     87double PolyOneShell(double dp[], double q); 
     88double PolyTwoShell(double dp[], double q); 
     89double PolyThreeShell(double dp[], double q); 
     90double PolyFourShell(double dp[], double q); 
     91double BCC_ParaCrystal(double dp[], double q); 
     92double FCC_ParaCrystal(double dp[], double q); 
     93double SC_ParaCrystal(double dp[], double q); 
    7494 
    7595//function prototypes 
    … …  
    106126double ThreeLevel(double dp[], double q); 
    107127double FourLevel(double dp[], double q); 
     128// 
     129double BroadPeak(double dp[], double q); 
     130double CorrLength(double dp[], double q); 
     131double TwoLorentzian(double dp[], double q); 
     132double TwoPowerLaw(double dp[], double q); 
     133double PolyGaussCoil(double dp[], double q); 
     134double GaussLorentzGel(double dp[], double q); 
     135double GaussianShell(double dp[], double q); 
     136 
    108137 
    109138// since the XOP and the library are separate chunks of compiled code 

  • sans/XOP_Dev/SANSAnalysis/lib/libSphere.c

    r235 r453  
    12241224} 
    12251225 
     1226double 
     1227OneShell(double dp[], double q) 
     1228{ 
     1229        // variables are: 
     1230        //[0] scale factor 
     1231        //[1] radius of core [] 
     1232        //[2] SLD of the core   [-2] 
     1233        //[3] thickness of the shell    [] 
     1234        //[4] SLD of the shell 
     1235        //[5] SLD of the solvent 
     1236        //[6] background        [cm-1] 
     1237         
     1238        double x,pi; 
     1239        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg;           //my local names 
     1240        double bes,f,vol,qr,contr,f2; 
     1241         
     1242        pi = 4.0*atan(1.0); 
     1243        x=q; 
     1244         
     1245        scale = dp[0]; 
     1246        rcore = dp[1]; 
     1247        rhocore = dp[2]; 
     1248        thick = dp[3]; 
     1249        rhoshel = dp[4]; 
     1250        rhosolv = dp[5]; 
     1251        bkg = dp[6]; 
     1252         
     1253        // core first, then add in shell 
     1254        qr=x*rcore; 
     1255        contr = rhocore-rhoshel; 
     1256        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1257        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
     1258        f = vol*bes*contr; 
     1259        //now the shell 
     1260        qr=x*(rcore+thick); 
     1261        contr = rhoshel-rhosolv; 
     1262        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1263        vol = 4.0*pi/3.0*pow((rcore+thick),3); 
     1264        f += vol*bes*contr; 
     1265         
     1266        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1267        f2 = f*f/vol*1.0e8; 
     1268         
     1269        //scale if desired 
     1270        f2 *= scale; 
     1271        // then add in the background 
     1272        f2 += bkg; 
     1273         
     1274        return(f2); 
     1275} 
     1276 
     1277double 
     1278TwoShell(double dp[], double q) 
     1279{ 
     1280        // variables are: 
     1281        //[0] scale factor 
     1282        //[1] radius of core [] 
     1283        //[2] SLD of the core   [-2] 
     1284        //[3] thickness of shell 1 [] 
     1285        //[4] SLD of shell 1 
     1286        //[5] thickness of shell 2 [] 
     1287        //[6] SLD of shell 2 
     1288        //[7] SLD of the solvent 
     1289        //[8] background        [cm-1] 
     1290         
     1291        double x,pi; 
     1292        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names 
     1293        double bes,f,vol,qr,contr,f2; 
     1294        double rhoshel2,thick2; 
     1295         
     1296        pi = 4.0*atan(1.0); 
     1297        x=q; 
     1298         
     1299        scale = dp[0]; 
     1300        rcore = dp[1]; 
     1301        rhocore = dp[2]; 
     1302        thick1 = dp[3]; 
     1303        rhoshel1 = dp[4]; 
     1304        thick2 = dp[5]; 
     1305        rhoshel2 = dp[6];        
     1306        rhosolv = dp[7]; 
     1307        bkg = dp[8]; 
     1308                // core first, then add in shells 
     1309        qr=x*rcore; 
     1310        contr = rhocore-rhoshel1; 
     1311        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1312        vol = 4.0*pi/3*rcore*rcore*rcore; 
     1313        f = vol*bes*contr; 
     1314        //now the shell (1) 
     1315        qr=x*(rcore+thick1); 
     1316        contr = rhoshel1-rhoshel2; 
     1317        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1318        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
     1319        f += vol*bes*contr; 
     1320        //now the shell (2) 
     1321        qr=x*(rcore+thick1+thick2); 
     1322        contr = rhoshel2-rhosolv; 
     1323        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1324        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
     1325        f += vol*bes*contr; 
     1326         
     1327                 
     1328        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1329        f2 = f*f/vol*1.0e8; 
     1330         
     1331        //scale if desired 
     1332        f2 *= scale; 
     1333        // then add in the background 
     1334        f2 += bkg; 
     1335         
     1336        return(f2); 
     1337} 
     1338 
     1339double 
     1340ThreeShell(double dp[], double q) 
     1341{ 
     1342        // variables are: 
     1343        //[0] scale factor 
     1344        //[1] radius of core [] 
     1345        //[2] SLD of the core   [-2] 
     1346        //[3] thickness of shell 1 [] 
     1347        //[4] SLD of shell 1 
     1348        //[5] thickness of shell 2 [] 
     1349        //[6] SLD of shell 2 
     1350        //[7] thickness of shell 3 
     1351        //[8] SLD of shell 3 
     1352        //[9] SLD of solvent 
     1353        //[10] background       [cm-1] 
     1354         
     1355        double x,pi; 
     1356        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names 
     1357        double bes,f,vol,qr,contr,f2; 
     1358        double rhoshel2,thick2,rhoshel3,thick3; 
     1359         
     1360        pi = 4.0*atan(1.0); 
     1361        x=q; 
     1362         
     1363        scale = dp[0]; 
     1364        rcore = dp[1]; 
     1365        rhocore = dp[2]; 
     1366        thick1 = dp[3]; 
     1367        rhoshel1 = dp[4]; 
     1368        thick2 = dp[5]; 
     1369        rhoshel2 = dp[6];        
     1370        thick3 = dp[7]; 
     1371        rhoshel3 = dp[8];        
     1372        rhosolv = dp[9]; 
     1373        bkg = dp[10]; 
     1374         
     1375                // core first, then add in shells 
     1376        qr=x*rcore; 
     1377        contr = rhocore-rhoshel1; 
     1378        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1379        vol = 4.0*pi/3*rcore*rcore*rcore; 
     1380        f = vol*bes*contr; 
     1381        //now the shell (1) 
     1382        qr=x*(rcore+thick1); 
     1383        contr = rhoshel1-rhoshel2; 
     1384        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1385        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
     1386        f += vol*bes*contr; 
     1387        //now the shell (2) 
     1388        qr=x*(rcore+thick1+thick2); 
     1389        contr = rhoshel2-rhoshel3; 
     1390        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1391        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
     1392        f += vol*bes*contr; 
     1393        //now the shell (3) 
     1394        qr=x*(rcore+thick1+thick2+thick3); 
     1395        contr = rhoshel3-rhosolv; 
     1396        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1397        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
     1398        f += vol*bes*contr; 
     1399                 
     1400        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1401        f2 = f*f/vol*1.0e8; 
     1402         
     1403        //scale if desired 
     1404        f2 *= scale; 
     1405        // then add in the background 
     1406        f2 += bkg; 
     1407         
     1408        return(f2); 
     1409} 
     1410 
     1411double 
     1412FourShell(double dp[], double q) 
     1413{ 
     1414        // variables are: 
     1415        //[0] scale factor 
     1416        //[1] radius of core [] 
     1417        //[2] SLD of the core   [-2] 
     1418        //[3] thickness of shell 1 [] 
     1419        //[4] SLD of shell 1 
     1420        //[5] thickness of shell 2 [] 
     1421        //[6] SLD of shell 2 
     1422        //[7] thickness of shell 3 
     1423        //[8] SLD of shell 3 
     1424        //[9] thickness of shell 3 
     1425        //[10] SLD of shell 3 
     1426        //[11] SLD of solvent 
     1427        //[12] background       [cm-1] 
     1428         
     1429        double x,pi; 
     1430        double scale,rcore,thick1,rhocore,rhoshel1,rhosolv,bkg;         //my local names 
     1431        double bes,f,vol,qr,contr,f2; 
     1432        double rhoshel2,thick2,rhoshel3,thick3,rhoshel4,thick4; 
     1433         
     1434        pi = 4.0*atan(1.0); 
     1435        x=q; 
     1436         
     1437        scale = dp[0]; 
     1438        rcore = dp[1]; 
     1439        rhocore = dp[2]; 
     1440        thick1 = dp[3]; 
     1441        rhoshel1 = dp[4]; 
     1442        thick2 = dp[5]; 
     1443        rhoshel2 = dp[6];        
     1444        thick3 = dp[7]; 
     1445        rhoshel3 = dp[8]; 
     1446        thick4 = dp[9]; 
     1447        rhoshel4 = dp[10];       
     1448        rhosolv = dp[11]; 
     1449        bkg = dp[12]; 
     1450         
     1451                // core first, then add in shells 
     1452        qr=x*rcore; 
     1453        contr = rhocore-rhoshel1; 
     1454        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1455        vol = 4.0*pi/3*rcore*rcore*rcore; 
     1456        f = vol*bes*contr; 
     1457        //now the shell (1) 
     1458        qr=x*(rcore+thick1); 
     1459        contr = rhoshel1-rhoshel2; 
     1460        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1461        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
     1462        f += vol*bes*contr; 
     1463        //now the shell (2) 
     1464        qr=x*(rcore+thick1+thick2); 
     1465        contr = rhoshel2-rhoshel3; 
     1466        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1467        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
     1468        f += vol*bes*contr; 
     1469        //now the shell (3) 
     1470        qr=x*(rcore+thick1+thick2+thick3); 
     1471        contr = rhoshel3-rhoshel4; 
     1472        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1473        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
     1474        f += vol*bes*contr; 
     1475        //now the shell (4) 
     1476        qr=x*(rcore+thick1+thick2+thick3+thick4); 
     1477        contr = rhoshel4-rhosolv; 
     1478        bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1479        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); 
     1480        f += vol*bes*contr; 
     1481         
     1482                 
     1483        // normalize to particle volume and rescale from [-1] to [cm-1] 
     1484        f2 = f*f/vol*1.0e8; 
     1485         
     1486        //scale if desired 
     1487        f2 *= scale; 
     1488        // then add in the background 
     1489        f2 += bkg; 
     1490         
     1491        return(f2); 
     1492} 
     1493 
     1494double 
     1495PolyOneShell(double dp[], double x) 
     1496{ 
     1497        double scale,rcore,thick,rhocore,rhoshel,rhosolv,bkg,pd,zz;             //my local names 
     1498        double va,vb,summ,yyy,zi; 
     1499        double answer,zp1,zp2,zp3,vpoly,range,temp_1sf[7],pi; 
     1500        int nord=76,ii; 
     1501         
     1502        pi = 4.0*atan(1.0); 
     1503         
     1504        scale = dp[0]; 
     1505        rcore = dp[1]; 
     1506        pd = dp[2]; 
     1507        rhocore = dp[3]; 
     1508        thick = dp[4]; 
     1509        rhoshel = dp[5]; 
     1510        rhosolv = dp[6]; 
     1511        bkg = dp[7]; 
     1512                 
     1513        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1514         
     1515        range = 8.0;                    //std deviations for the integration 
     1516        va = rcore*(1.0-range*pd); 
     1517        if (va<0) { 
     1518                va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1519        } 
     1520        if (pd>0.3) { 
     1521                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1522        } 
     1523        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1524         
     1525        //temp set scale=1 and bkg=0 for quadrature calc 
     1526        temp_1sf[0] = 1.0; 
     1527        temp_1sf[1] = dp[1];    //the core radius will be changed in the loop 
     1528        temp_1sf[2] = dp[3]; 
     1529        temp_1sf[3] = dp[4]; 
     1530        temp_1sf[4] = dp[5]; 
     1531        temp_1sf[5] = dp[6]; 
     1532        temp_1sf[6] = 0.0; 
     1533         
     1534        summ = 0.0;             // initialize integral 
     1535        for(ii=0;ii<nord;ii+=1) { 
     1536                // calculate Gauss points on integration interval (r-value for evaluation) 
     1537                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1538                temp_1sf[1] = zi; 
     1539                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * OneShell(temp_1sf,x); 
     1540                //un-normalize by volume 
     1541                yyy *= 4.0*pi/3.0*pow((zi+thick),3); 
     1542                summ += yyy;            //add to the running total of the quadrature 
     1543        } 
     1544        // calculate value of integral to return 
     1545        answer = (vb-va)/2.0*summ; 
     1546         
     1547        //re-normalize by the average volume 
     1548        zp1 = zz + 1.0; 
     1549        zp2 = zz + 2.0; 
     1550        zp3 = zz + 3.0; 
     1551        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick),3); 
     1552        answer /= vpoly; 
     1553//scale 
     1554        answer *= scale; 
     1555// add in the background 
     1556        answer += bkg; 
     1557                 
     1558    return(answer); 
     1559} 
     1560 
     1561double 
     1562PolyTwoShell(double dp[], double x) 
     1563{ 
     1564        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names 
     1565        double va,vb,summ,yyy,zi; 
     1566        double answer,zp1,zp2,zp3,vpoly,range,temp_2sf[9],pi; 
     1567        int nord=76,ii; 
     1568        double thick1,thick2; 
     1569        double rhoshel1,rhoshel2; 
     1570         
     1571        scale = dp[0]; 
     1572        rcore = dp[1]; 
     1573        pd = dp[2]; 
     1574        rhocore = dp[3]; 
     1575        thick1 = dp[4]; 
     1576        rhoshel1 = dp[5]; 
     1577        thick2 = dp[6]; 
     1578        rhoshel2 = dp[7]; 
     1579        rhosolv = dp[8]; 
     1580        bkg = dp[9]; 
     1581         
     1582        pi = 4.0*atan(1.0); 
     1583                 
     1584        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1585         
     1586        range = 8.0;                    //std deviations for the integration 
     1587        va = rcore*(1.0-range*pd); 
     1588        if (va<0) { 
     1589                va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1590        } 
     1591        if (pd>0.3) { 
     1592                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1593        } 
     1594        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1595         
     1596        //temp set scale=1 and bkg=0 for quadrature calc 
     1597        temp_2sf[0] = 1.0; 
     1598        temp_2sf[1] = dp[1];            //the core radius will be changed in the loop 
     1599        temp_2sf[2] = dp[3]; 
     1600        temp_2sf[3] = dp[4]; 
     1601        temp_2sf[4] = dp[5]; 
     1602        temp_2sf[5] = dp[6]; 
     1603        temp_2sf[6] = dp[7]; 
     1604        temp_2sf[7] = dp[8]; 
     1605        temp_2sf[8] = 0.0; 
     1606         
     1607        summ = 0.0;             // initialize integral 
     1608        for(ii=0;ii<nord;ii+=1) { 
     1609                // calculate Gauss points on integration interval (r-value for evaluation) 
     1610                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1611                temp_2sf[1] = zi; 
     1612                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * TwoShell(temp_2sf,x); 
     1613                //un-normalize by volume 
     1614                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2),3); 
     1615                summ += yyy;            //add to the running total of the quadrature 
     1616        } 
     1617        // calculate value of integral to return 
     1618        answer = (vb-va)/2.0*summ; 
     1619         
     1620        //re-normalize by the average volume 
     1621        zp1 = zz + 1.0; 
     1622        zp2 = zz + 2.0; 
     1623        zp3 = zz + 3.0; 
     1624        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2),3); 
     1625        answer /= vpoly; 
     1626//scale 
     1627        answer *= scale; 
     1628// add in the background 
     1629        answer += bkg; 
     1630                 
     1631    return(answer); 
     1632} 
     1633 
     1634double 
     1635PolyThreeShell(double dp[], double x) 
     1636{ 
     1637        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names 
     1638        double va,vb,summ,yyy,zi; 
     1639        double answer,zp1,zp2,zp3,vpoly,range,temp_3sf[11],pi; 
     1640        int nord=76,ii; 
     1641        double thick1,thick2,thick3; 
     1642        double rhoshel1,rhoshel2,rhoshel3; 
     1643         
     1644        scale = dp[0]; 
     1645        rcore = dp[1]; 
     1646        pd = dp[2]; 
     1647        rhocore = dp[3]; 
     1648        thick1 = dp[4]; 
     1649        rhoshel1 = dp[5]; 
     1650        thick2 = dp[6]; 
     1651        rhoshel2 = dp[7]; 
     1652        thick3 = dp[8]; 
     1653        rhoshel3 = dp[9]; 
     1654        rhosolv = dp[10]; 
     1655        bkg = dp[11]; 
     1656         
     1657        pi = 4.0*atan(1.0); 
     1658                 
     1659        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1660         
     1661        range = 8.0;                    //std deviations for the integration 
     1662        va = rcore*(1.0-range*pd); 
     1663        if (va<0) { 
     1664                va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1665        } 
     1666        if (pd>0.3) { 
     1667                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1668        } 
     1669        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1670         
     1671        //temp set scale=1 and bkg=0 for quadrature calc 
     1672        temp_3sf[0] = 1.0; 
     1673        temp_3sf[1] = dp[1];            //the core radius will be changed in the loop 
     1674        temp_3sf[2] = dp[3]; 
     1675        temp_3sf[3] = dp[4]; 
     1676        temp_3sf[4] = dp[5]; 
     1677        temp_3sf[5] = dp[6]; 
     1678        temp_3sf[6] = dp[7]; 
     1679        temp_3sf[7] = dp[8]; 
     1680        temp_3sf[8] = dp[9]; 
     1681        temp_3sf[9] = dp[10]; 
     1682        temp_3sf[10] = 0.0; 
     1683         
     1684        summ = 0.0;             // initialize integral 
     1685        for(ii=0;ii<nord;ii+=1) { 
     1686                // calculate Gauss points on integration interval (r-value for evaluation) 
     1687                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1688                temp_3sf[1] = zi; 
     1689                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * ThreeShell(temp_3sf,x); 
     1690                //un-normalize by volume 
     1691                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3),3); 
     1692                summ += yyy;            //add to the running total of the quadrature 
     1693        } 
     1694        // calculate value of integral to return 
     1695        answer = (vb-va)/2.0*summ; 
     1696         
     1697        //re-normalize by the average volume 
     1698        zp1 = zz + 1.0; 
     1699        zp2 = zz + 2.0; 
     1700        zp3 = zz + 3.0; 
     1701        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3),3); 
     1702        answer /= vpoly; 
     1703//scale 
     1704        answer *= scale; 
     1705// add in the background 
     1706        answer += bkg; 
     1707                 
     1708    return(answer); 
     1709} 
     1710 
     1711double 
     1712PolyFourShell(double dp[], double x) 
     1713{ 
     1714        double scale,rcore,rhocore,rhosolv,bkg,pd,zz;           //my local names 
     1715        double va,vb,summ,yyy,zi; 
     1716        double answer,zp1,zp2,zp3,vpoly,range,temp_4sf[13],pi; 
     1717        int nord=76,ii; 
     1718        double thick1,thick2,thick3,thick4; 
     1719        double rhoshel1,rhoshel2,rhoshel3,rhoshel4; 
     1720         
     1721        scale = dp[0]; 
     1722        rcore = dp[1]; 
     1723        pd = dp[2]; 
     1724        rhocore = dp[3]; 
     1725        thick1 = dp[4]; 
     1726        rhoshel1 = dp[5]; 
     1727        thick2 = dp[6]; 
     1728        rhoshel2 = dp[7]; 
     1729        thick3 = dp[8]; 
     1730        rhoshel3 = dp[9]; 
     1731        thick4 = dp[10]; 
     1732        rhoshel4 = dp[11]; 
     1733        rhosolv = dp[12]; 
     1734        bkg = dp[13]; 
     1735         
     1736        pi = 4.0*atan(1.0); 
     1737                 
     1738        zz = (1.0/pd)*(1.0/pd)-1.0;             //polydispersity of the core only 
     1739         
     1740        range = 8.0;                    //std deviations for the integration 
     1741        va = rcore*(1.0-range*pd); 
     1742        if (va<0) { 
     1743                va=0;           //otherwise numerical error when pd >= 0.3, making a<0 
     1744        } 
     1745        if (pd>0.3) { 
     1746                range = range + (pd-0.3)*18.0;          //stretch upper range to account for skewed tail 
     1747        } 
     1748        vb = rcore*(1.0+range*pd);              // is this far enough past avg radius? 
     1749         
     1750        //temp set scale=1 and bkg=0 for quadrature calc 
     1751        temp_4sf[0] = 1.0; 
     1752        temp_4sf[1] = dp[1];            //the core radius will be changed in the loop 
     1753        temp_4sf[2] = dp[3]; 
     1754        temp_4sf[3] = dp[4]; 
     1755        temp_4sf[4] = dp[5]; 
     1756        temp_4sf[5] = dp[6]; 
     1757        temp_4sf[6] = dp[7]; 
     1758        temp_4sf[7] = dp[8]; 
     1759        temp_4sf[8] = dp[9]; 
     1760        temp_4sf[9] = dp[10]; 
     1761        temp_4sf[10] = dp[11]; 
     1762        temp_4sf[11] = dp[12]; 
     1763        temp_4sf[12] = 0.0; 
     1764                 
     1765        summ = 0.0;             // initialize integral 
     1766        for(ii=0;ii<nord;ii+=1) { 
     1767                // calculate Gauss points on integration interval (r-value for evaluation) 
     1768                zi = ( Gauss76Z[ii]*(vb-va) + vb + va )/2.0; 
     1769                temp_4sf[1] = zi; 
     1770                yyy = Gauss76Wt[ii] * SchulzPoint(zi,rcore,zz) * FourShell(temp_4sf,x); 
     1771                //un-normalize by volume 
     1772                yyy *= 4.0*pi/3.0*pow((zi+thick1+thick2+thick3+thick4),3); 
     1773                summ += yyy;            //add to the running total of the quadrature 
     1774        } 
     1775        // calculate value of integral to return 
     1776        answer = (vb-va)/2.0*summ; 
     1777         
     1778        //re-normalize by the average volume 
     1779        zp1 = zz + 1.0; 
     1780        zp2 = zz + 2.0; 
     1781        zp3 = zz + 3.0; 
     1782        vpoly = 4.0*pi/3.0*zp3*zp2/zp1/zp1*pow((rcore+thick1+thick2+thick3+thick4),3); 
     1783        answer /= vpoly; 
     1784//scale 
     1785        answer *= scale; 
     1786// add in the background 
     1787        answer += bkg; 
     1788                 
     1789    return(answer); 
     1790} 
     1791 
     1792 
     1793/*      BCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x 
     1794 
     1795Uses 150 pt Gaussian quadrature for both integrals 
     1796 
     1797*/ 
     1798double 
     1799BCC_ParaCrystal(double w[], double x) 
     1800{ 
     1801        int i,j; 
     1802        double Pi; 
     1803        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave 
     1804        int nordi=150;                  //order of integration 
     1805        int nordj=150; 
     1806        double va,vb;           //upper and lower integration limits 
     1807        double summ,zi,yyy,answer;                      //running tally of integration 
     1808        double summj,vaj,vbj,zij;                       //for the inner integration 
     1809         
     1810        Pi = 4.0*atan(1.0); 
     1811        va = 0.0; 
     1812        vb = 2.0*Pi;            //orintational average, outer integral 
     1813        vaj = 0.0; 
     1814        vbj = Pi;               //endpoints of inner integral 
     1815         
     1816        summ = 0.0;                     //initialize intergral 
     1817         
     1818        scale = w[0]; 
     1819        Dnn = w[1];                                     //Nearest neighbor distance A 
     1820        gg = w[2];                                      //Paracrystal distortion factor 
     1821        Rad = w[3];                                     //Sphere radius 
     1822        sld = w[4]; 
     1823        sldSolv = w[5]; 
     1824        background = w[6];  
     1825         
     1826        contrast = sld - sldSolv; 
     1827         
     1828        //Volume fraction calculated from lattice symmetry and sphere radius 
     1829        latticeScale = 2.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn/sqrt(3.0/4.0),3); 
     1830         
     1831        for(i=0;i<nordi;i++) { 
     1832                //setup inner integral over the ellipsoidal cross-section 
     1833                summj=0.0; 
     1834                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi 
     1835                for(j=0;j<nordj;j++) { 
     1836                        //20 gauss points for the inner integral 
     1837                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta 
     1838                        yyy = Gauss150Wt[j] * BCC_Integrand(w,x,zi,zij); 
     1839                        summj += yyy; 
     1840                } 
     1841                //now calculate the value of the inner integral 
     1842                answer = (vbj-vaj)/2.0*summj; 
     1843                 
     1844                //now calculate outer integral 
     1845                yyy = Gauss150Wt[i] * answer; 
     1846                summ += yyy; 
     1847        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     1848         
     1849        answer = (vb-va)/2.0*summ; 
     1850        // Multiply by contrast^2 
     1851        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 
     1852        // add in the background 
     1853        answer += background; 
     1854         
     1855        return answer; 
     1856} 
     1857 
     1858// xx is phi (outer) 
     1859// yy is theta (inner) 
     1860double 
     1861BCC_Integrand(double w[], double qq, double xx, double yy) { 
     1862         
     1863        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi; 
     1864         
     1865        Dnn = w[1]; //Nearest neighbor distance A 
     1866        gg = w[2]; //Paracrystal distortion factor 
     1867        aa = Dnn; 
     1868        Da = gg*aa; 
     1869         
     1870        Pi = 4.0*atan(1.0); 
     1871        temp1 = qq*qq*Da*Da; 
     1872        temp3 = qq*aa;   
     1873         
     1874        retVal = BCCeval(yy,xx,temp1,temp3); 
     1875        retVal /=4.0*Pi; 
     1876         
     1877        return(retVal); 
     1878} 
     1879 
     1880double 
     1881BCCeval(double Theta, double Phi, double temp1, double temp3) { 
     1882 
     1883        double temp6,temp7,temp8,temp9,temp10; 
     1884        double result; 
     1885         
     1886        temp6 = sin(Theta); 
     1887        temp7 = sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)+cos(Theta); 
     1888        temp8 = -1.0*sin(Theta)*cos(Phi)-sin(Theta)*sin(Phi)+cos(Theta); 
     1889        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi)-cos(Theta); 
     1890        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
     1891        result = pow(1.0-(temp10*temp10),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 
     1892         
     1893        return (result); 
     1894} 
     1895 
     1896double 
     1897SphereForm_Paracrystal(double radius, double delrho, double x) {                                         
     1898         
     1899        double bes,f,vol,f2,pi; 
     1900        pi = 4.0*atan(1.0); 
     1901        // 
     1902        //handle q==0 separately 
     1903        if(x==0) { 
     1904                f = 4.0/3.0*pi*radius*radius*radius*delrho*delrho*1.0e8; 
     1905                return(f); 
     1906        } 
     1907         
     1908        bes = 3.0*(sin(x*radius)-x*radius*cos(x*radius))/(x*x*x)/(radius*radius*radius); 
     1909        vol = 4.0*pi/3.0*radius*radius*radius; 
     1910        f = vol*bes*delrho      ;       // [=]  
     1911        // normalize to single particle volume, convert to 1/cm 
     1912        f2 = f * f / vol * 1.0e8;               // [=] 1/cm 
     1913         
     1914        return (f2); 
     1915} 
     1916 
     1917/*      FCC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x 
     1918 
     1919Uses 150 pt Gaussian quadrature for both integrals 
     1920 
     1921*/ 
     1922double 
     1923FCC_ParaCrystal(double w[], double x) 
     1924{ 
     1925        int i,j; 
     1926        double Pi; 
     1927        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave 
     1928        int nordi=150;                  //order of integration 
     1929        int nordj=150; 
     1930        double va,vb;           //upper and lower integration limits 
     1931        double summ,zi,yyy,answer;                      //running tally of integration 
     1932        double summj,vaj,vbj,zij;                       //for the inner integration 
     1933         
     1934        Pi = 4.0*atan(1.0); 
     1935        va = 0.0; 
     1936        vb = 2.0*Pi;            //orintational average, outer integral 
     1937        vaj = 0.0; 
     1938        vbj = Pi;               //endpoints of inner integral 
     1939         
     1940        summ = 0.0;                     //initialize intergral 
     1941         
     1942        scale = w[0]; 
     1943        Dnn = w[1];                                     //Nearest neighbor distance A 
     1944        gg = w[2];                                      //Paracrystal distortion factor 
     1945        Rad = w[3];                                     //Sphere radius 
     1946        sld = w[4]; 
     1947        sldSolv = w[5]; 
     1948        background = w[6];  
     1949         
     1950        contrast = sld - sldSolv; 
     1951        //Volume fraction calculated from lattice symmetry and sphere radius 
     1952        latticeScale = 4.0*(4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn*sqrt(2.0),3); 
     1953         
     1954        for(i=0;i<nordi;i++) { 
     1955                //setup inner integral over the ellipsoidal cross-section 
     1956                summj=0.0; 
     1957                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi 
     1958                for(j=0;j<nordj;j++) { 
     1959                        //20 gauss points for the inner integral 
     1960                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta 
     1961                        yyy = Gauss150Wt[j] * FCC_Integrand(w,x,zi,zij); 
     1962                        summj += yyy; 
     1963                } 
     1964                //now calculate the value of the inner integral 
     1965                answer = (vbj-vaj)/2.0*summj; 
     1966                 
     1967                //now calculate outer integral 
     1968                yyy = Gauss150Wt[i] * answer; 
     1969                summ += yyy; 
     1970        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     1971         
     1972        answer = (vb-va)/2.0*summ; 
     1973        // Multiply by contrast^2 
     1974        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 
     1975        // add in the background 
     1976        answer += background; 
     1977         
     1978        return answer; 
     1979} 
     1980 
     1981 
     1982// xx is phi (outer) 
     1983// yy is theta (inner) 
     1984double 
     1985FCC_Integrand(double w[], double qq, double xx, double yy) { 
     1986         
     1987        double retVal,temp1,temp3,aa,Da,Dnn,gg,Pi; 
     1988         
     1989        Pi = 4.0*atan(1.0); 
     1990        Dnn = w[1]; //Nearest neighbor distance A 
     1991        gg = w[2]; //Paracrystal distortion factor 
     1992        aa = Dnn; 
     1993        Da = gg*aa; 
     1994         
     1995        temp1 = qq*qq*Da*Da; 
     1996        temp3 = qq*aa; 
     1997         
     1998        retVal = FCCeval(yy,xx,temp1,temp3); 
     1999        retVal /=4*Pi; 
     2000         
     2001        return(retVal); 
     2002} 
     2003 
     2004double 
     2005FCCeval(double Theta, double Phi, double temp1, double temp3) { 
     2006 
     2007        double temp6,temp7,temp8,temp9,temp10; 
     2008        double result; 
     2009         
     2010        temp6 = sin(Theta); 
     2011        temp7 = sin(Theta)*sin(Phi)+cos(Theta); 
     2012        temp8 = -1.0*sin(Theta)*cos(Phi)+cos(Theta); 
     2013        temp9 = -1.0*sin(Theta)*cos(Phi)+sin(Theta)*sin(Phi); 
     2014        temp10 = exp((-1.0/8.0)*temp1*((temp7*temp7)+(temp8*temp8)+(temp9*temp9))); 
     2015        result = pow((1.0-(temp10*temp10)),3)*temp6/((1.0-2.0*temp10*cos(0.5*temp3*(temp7))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp8))+(temp10*temp10))*(1.0-2.0*temp10*cos(0.5*temp3*(temp9))+(temp10*temp10))); 
     2016         
     2017        return (result); 
     2018} 
     2019 
     2020 
     2021/*      SC_ParaCrystal  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x 
     2022 
     2023Uses 150 pt Gaussian quadrature for both integrals 
     2024 
     2025*/ 
     2026double 
     2027SC_ParaCrystal(double w[], double x) 
     2028{ 
     2029        int i,j; 
     2030        double Pi; 
     2031        double scale,Dnn,gg,Rad,contrast,background,latticeScale,sld,sldSolv;           //local variables of coefficient wave 
     2032        int nordi=150;                  //order of integration 
     2033        int nordj=150; 
     2034        double va,vb;           //upper and lower integration limits 
     2035        double summ,zi,yyy,answer;                      //running tally of integration 
     2036        double summj,vaj,vbj,zij;                       //for the inner integration 
     2037         
     2038        Pi = 4.0*atan(1.0); 
     2039        va = 0.0; 
     2040        vb = 2.0*Pi;            //orintational average, outer integral 
     2041        vaj = 0.0; 
     2042        vbj = Pi;               //endpoints of inner integral 
     2043         
     2044        summ = 0.0;                     //initialize intergral 
     2045         
     2046        scale = w[0]; 
     2047        Dnn = w[1];                                     //Nearest neighbor distance A 
     2048        gg = w[2];                                      //Paracrystal distortion factor 
     2049        Rad = w[3];                                     //Sphere radius 
     2050        sld = w[4]; 
     2051        sldSolv = w[5]; 
     2052        background = w[6];  
     2053         
     2054        contrast = sld - sldSolv; 
     2055        //Volume fraction calculated from lattice symmetry and sphere radius 
     2056        latticeScale = (4.0/3.0)*Pi*(Rad*Rad*Rad)/pow(Dnn,3); 
     2057         
     2058        for(i=0;i<nordi;i++) { 
     2059                //setup inner integral over the ellipsoidal cross-section 
     2060                summj=0.0; 
     2061                zi = ( Gauss150Z[i]*(vb-va) + va + vb )/2.0;            //the outer dummy is phi 
     2062                for(j=0;j<nordj;j++) { 
     2063                        //20 gauss points for the inner integral 
     2064                        zij = ( Gauss150Z[j]*(vbj-vaj) + vaj + vbj )/2.0;               //the inner dummy is theta 
     2065                        yyy = Gauss150Wt[j] * SC_Integrand(w,x,zi,zij); 
     2066                        summj += yyy; 
     2067                } 
     2068                //now calculate the value of the inner integral 
     2069                answer = (vbj-vaj)/2.0*summj; 
     2070                 
     2071                //now calculate outer integral 
     2072                yyy = Gauss150Wt[i] * answer; 
     2073                summ += yyy; 
     2074        }               //final scaling is done at the end of the function, after the NT_FP64 case 
     2075         
     2076        answer = (vb-va)/2.0*summ; 
     2077        // Multiply by contrast^2 
     2078        answer *= SphereForm_Paracrystal(Rad,contrast,x)*scale*latticeScale; 
     2079        // add in the background 
     2080        answer += background; 
     2081         
     2082        return answer; 
     2083} 
     2084 
     2085// xx is phi (outer) 
     2086// yy is theta (inner) 
     2087double 
     2088SC_Integrand(double w[], double qq, double xx, double yy) { 
     2089         
     2090        double retVal,temp1,temp2,temp3,temp4,temp5,aa,Da,Dnn,gg,Pi; 
     2091         
     2092        Pi = 4.0*atan(1.0); 
     2093        Dnn = w[1]; //Nearest neighbor distance A 
     2094        gg = w[2]; //Paracrystal distortion factor 
     2095        aa = Dnn; 
     2096        Da = gg*aa; 
     2097         
     2098        temp1 = qq*qq*Da*Da; 
     2099        temp2 = pow( 1.0-exp(-1.0*temp1) ,3); 
     2100        temp3 = qq*aa; 
     2101        temp4 = 2.0*exp(-0.5*temp1); 
     2102        temp5 = exp(-1.0*temp1); 
     2103         
     2104         
     2105        retVal = temp2*SCeval(yy,xx,temp3,temp4,temp5); 
     2106        retVal /= 4*Pi; 
     2107         
     2108        return(retVal); 
     2109} 
     2110 
     2111double 
     2112SCeval(double Theta, double Phi, double temp3, double temp4, double temp5) { //Function to calculate integrand values for simple cubic structure 
     2113 
     2114        double temp6,temp7,temp8,temp9; //Theta and phi dependent parts of the equation 
     2115        double result; 
     2116         
     2117        temp6 = sin(Theta); 
     2118        temp7 = -1.0*temp3*sin(Theta)*cos(Phi); 
     2119        temp8 = temp3*sin(Theta)*sin(Phi); 
     2120        temp9 = temp3*cos(Theta); 
     2121        result = temp6/((1.0-temp4*cos((temp7))+temp5)*(1.0-temp4*cos((temp8))+temp5)*(1.0-temp4*cos((temp9))+temp5)); 
     2122         
     2123        return (result); 
     2124} 
     2125 

  • sans/XOP_Dev/SANSAnalysis/lib/libSphere.h

    r101 r453  
    1919double BinaryHS_PSF12(double dp[], double q); 
    2020double BinaryHS_PSF22(double dp[], double q); 
     21double OneShell(double dp[], double q); 
     22double TwoShell(double dp[], double q); 
     23double ThreeShell(double dp[], double q); 
     24double FourShell(double dp[], double q); 
     25double PolyOneShell(double dp[], double q); 
     26double PolyTwoShell(double dp[], double q); 
     27double PolyThreeShell(double dp[], double q); 
     28double PolyFourShell(double dp[], double q); 
     29double BCC_ParaCrystal(double dp[], double q); 
     30double FCC_ParaCrystal(double dp[], double q); 
     31double SC_ParaCrystal(double dp[], double q); 
    2132 
    2233//function prototypes 
    … …  
    2738double SchulzSphere_Fn(double scale, double ravg, double pd, double rho, double rhos, double x); 
    2839int ashcroft(double qval, double r2, double nf2, double aa, double phi, double *s11, double *s22, double *s12); 
     40double BCC_Integrand(double w[], double qq, double xx, double yy); 
     41double BCCeval(double Theta, double Phi, double temp1, double temp3); 
     42double SphereForm_Paracrystal(double radius, double delrho, double x); 
     43double FCC_Integrand(double w[], double qq, double xx, double yy); 
     44double FCCeval(double Theta, double Phi, double temp1, double temp3); 
     45double SC_Integrand(double w[], double qq, double xx, double yy); 
     46double SCeval(double Theta, double Phi, double temp3, double temp4, double temp5); 
     47 
     48 
    2949 
    3050static double SchulzPoint(double x, double avg, double zz); 
    … …  
    3252static double Gauss_distr(double sig, double avg, double pt); 
    3353static double LogNormal_distr(double sig, double mu, double pt); 
     54 

  • sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.c

    r97 r453  
    285285} 
    286286 
     287double 
     288BroadPeak(double dp[], double q) 
     289{ 
     290        // variables are:                                                        
     291        //[0] Porod term scaling 
     292        //[1] Porod exponent 
     293        //[2] Lorentzian term scaling 
     294        //[3] Lorentzian screening length [A] 
     295        //[4] peak location [1/A] 
     296        //[5] Lorentzian exponent 
     297        //[6] background 
     298         
     299        double aa,nn,cc,LL,Qzero,mm,bgd,inten,qval; 
     300        qval= q; 
     301        aa = dp[0]; 
     302        nn = dp[1]; 
     303        cc = dp[2];  
     304        LL = dp[3];  
     305        Qzero = dp[4];  
     306        mm = dp[5];  
     307        bgd = dp[6];  
     308         
     309        inten = aa/pow(qval,nn); 
     310        inten += cc/(1.0 + pow((fabs(qval-Qzero)*LL),mm) ); 
     311        inten += bgd; 
     312         
     313        return(inten); 
     314} 
     315 
     316double 
     317CorrLength(double dp[], double q) 
     318{ 
     319        // variables are:                                                        
     320        //[0] Porod term scaling 
     321        //[1] Porod exponent 
     322        //[2] Lorentzian term scaling 
     323        //[3] Lorentzian screening length [A] 
     324        //[4] Lorentzian exponent 
     325        //[5] background 
     326         
     327        double aa,nn,cc,LL,mm,bgd,inten,qval; 
     328        qval= q; 
     329        aa = dp[0]; 
     330        nn = dp[1]; 
     331        cc = dp[2];  
     332        LL = dp[3];  
     333        mm = dp[4];  
     334        bgd = dp[5];  
     335         
     336        inten = aa/pow(qval,nn); 
     337        inten += cc/(1.0 + pow((qval*LL),mm) ); 
     338        inten += bgd; 
     339         
     340        return(inten); 
     341} 
     342 
     343double 
     344TwoLorentzian(double dp[], double q) 
     345{ 
     346        // variables are:                                                        
     347        //[0] Lorentzian term scaling 
     348        //[1] Lorentzian screening length [A] 
     349        //[2] Lorentzian exponent 
     350        //[3] Lorentzian #2 term scaling 
     351        //[4] Lorentzian #2 screening length [A] 
     352        //[5] Lorentzian #2 exponent 
     353        //[6] background 
     354                 
     355        double aa,LL1,nn,cc,LL2,mm,bgd,inten,qval; 
     356        qval= q; 
     357        aa = dp[0]; 
     358        LL1 = dp[1]; 
     359        nn = dp[2];  
     360        cc = dp[3];  
     361        LL2 = dp[4];  
     362        mm = dp[5];  
     363        bgd= dp[6]; 
     364         
     365        inten = aa/(1.0 + pow((qval*LL1),nn) ); 
     366        inten += cc/(1.0 + pow((qval*LL2),mm) ); 
     367        inten += bgd; 
     368         
     369        return(inten); 
     370} 
     371 
     372double 
     373TwoPowerLaw(double dp[], double q) 
     374{ 
     375        //[0] Coefficient 
     376        //[1] (-) Power @ low Q 
     377        //[2] (-) Power @ high Q 
     378        //[3] crossover Q-value 
     379        //[4] incoherent background 
     380                 
     381        double A, m1,m2,qc,bgd,scale,inten,qval; 
     382        qval= q; 
     383        A = dp[0]; 
     384        m1 = dp[1]; 
     385        m2 = dp[2];  
     386        qc = dp[3];  
     387        bgd = dp[4];  
     388         
     389        if(qval<=qc){ 
     390                inten = A*pow(qval,-1.0*m1); 
     391        } else { 
     392                scale = A*pow(qc,-1.0*m1) / pow(qc,-1.0*m2); 
     393                inten = scale*pow(qval,-1.0*m2); 
     394        } 
     395         
     396        inten += bgd; 
     397         
     398        return(inten); 
     399} 
     400 
     401double 
     402PolyGaussCoil(double dp[], double x) 
     403{ 
     404        //w[0] = scale 
     405        //w[1] = radius of gyration [] 
     406        //w[2] = polydispersity, ratio of Mw/Mn 
     407        //w[3] = bkg [cm-1] 
     408                 
     409        double scale,bkg,Rg,uval,Mw_Mn,inten,xi; 
     410 
     411        scale = dp[0]; 
     412        Rg = dp[1]; 
     413        Mw_Mn = dp[2];  
     414        bkg = dp[3];  
     415         
     416        uval = Mw_Mn - 1.0; 
     417        if(uval == 0) { 
     418                uval = 1e-6;            //avoid divide by zero error 
     419        } 
     420         
     421        xi = Rg*Rg*x*x/(1.0+2.0*uval); 
     422         
     423        if(xi < 1e-3) { 
     424                return(scale+bkg);              //limiting value 
     425        } 
     426         
     427        inten = 2.0*(pow((1.0+uval*xi),(-1.0/uval))+xi-1.0); 
     428        inten /= (1.0+uval)*xi*xi; 
     429 
     430        inten *= scale; 
     431        //add in the background 
     432        inten += bkg;    
     433        return(inten); 
     434} 
     435 
     436double 
     437GaussLorentzGel(double dp[], double x) 
     438{ 
     439        //[0] Gaussian scale factor 
     440        //[1] Gaussian (static) screening length 
     441        //[2] Lorentzian (fluctuation) scale factor 
     442        //[3] Lorentzian screening length 
     443        //[4] incoherent background 
     444                 
     445        double Ig0,gg,Il0,ll,bgd,inten; 
     446 
     447        Ig0 = dp[0]; 
     448        gg = dp[1]; 
     449        Il0 = dp[2];  
     450        ll = dp[3]; 
     451        bgd = dp[4];  
     452         
     453        inten = Ig0*exp(-1.0*x*x*gg*gg/2.0) + Il0/(1 + (x*ll)*(x*ll)) + bgd; 
     454         
     455        return(inten); 
     456} 
     457 
    287458 
    288459static double 
    … …  
    303474} 
    304475 
    305  
     476double 
     477GaussianShell(double w[], double x) 
     478{ 
     479        // variables are:                                                        
     480        //[0] scale 
     481        //[1] radius () 
     482        //[2] thick () (thickness parameter - this is the std. dev. of the Gaussian width of the shell) 
     483        //[3] polydispersity of the radius 
     484        //[4] sld shell (-2) 
     485        //[5] sld solvent 
     486        //[6] background (cm-1) 
     487         
     488        double scale,rad,delrho,bkg,del,thick,pd,sig,pi; 
     489        double t1,t2,t3,t4,retval,exfact,vshell,vexcl,sldShell,sldSolvent; 
     490        scale = w[0]; 
     491        rad = w[1]; 
     492        thick = w[2]; 
     493        pd = w[3]; 
     494        sldShell = w[4]; 
     495        sldSolvent = w[5]; 
     496        bkg = w[6]; 
     497         
     498        delrho = w[4] - w[5]; 
     499        sig = pd*rad; 
     500         
     501        pi = 4.0*atan(1.0); 
     502         
     503        ///APPROXIMATION (see eqn 4 - but not a bad approximation) 
     504        // del is the equivalent shell thickness with sharp boundaries, centered at mean radius 
     505        del = thick*sqrt(2.0*pi); 
     506         
     507        // calculate the polydisperse shell volume and the excluded volume 
     508        vshell=4.0*pi/3.0*( pow((rad+del/2.0),3) - pow((rad-del/2.0),3) ) *(1.0+pd*pd); 
     509        vexcl=4.0*pi/3.0*( pow((rad+del/2.0),3) ) *(1.0+pd*pd); 
     510         
     511        //intensity, eqn 9(a-d) 
     512        exfact = exp(-2.0*sig*sig*x*x); 
     513         
     514        t1 = 0.5*x*x*thick*thick*thick*thick*(1.0+cos(2.0*x*rad)*exfact); 
     515        t2 = x*thick*thick*(rad*sin(2.0*x*rad) + 2.0*x*sig*sig*cos(2.0*x*rad))*exfact; 
     516        t3 = 0.5*rad*rad*(1.0-cos(2.0*x*rad)*exfact); 
     517        t4 = 0.5*sig*sig*(1.0+4.0*x*rad*sin(2.0*x*rad)*exfact+cos(2.0*x*rad)*(4.0*sig*sig*x*x-1.0)*exfact); 
     518         
     519        retval = t1+t2+t3+t4; 
     520        retval *= exp(-1.0*x*x*thick*thick); 
     521        retval *= (del*del/x/x); 
     522        retval *= 16.0*pi*pi*delrho*delrho*scale; 
     523        retval *= 1.0e8; 
     524         
     525        //NORMALIZED by the AVERAGE shell volume, since scale is the volume fraction of material 
     526//      retval /= vshell 
     527        retval /= vexcl; 
     528        //re-normalize by polydisperse sphere volume, Gaussian distribution 
     529        retval /= (1.0+3.0*pd*pd); 
     530         
     531        retval += bkg; 
     532         
     533        return(retval); 
     534} 
     535 
     536 

  • sans/XOP_Dev/SANSAnalysis/lib/libTwoPhase.h

    r97 r453  
    1414double ThreeLevel(double dp[], double q); 
    1515double FourLevel(double dp[], double q); 
     16double BroadPeak(double dp[], double q); 
     17double CorrLength(double dp[], double q); 
     18double TwoLorentzian(double dp[], double q); 
     19double TwoPowerLaw(double dp[], double q); 
     20double PolyGaussCoil(double dp[], double q); 
     21double GaussLorentzGel(double dp[], double q); 
     22double GaussianShell(double dp[], double q); 
     23 
    1624 
    1725/* internal functions */ 
Note: See TracChangeset for help on using the changeset viewer.