Changeset 201af9f in sasview


Ignore:
Timestamp:
Apr 5, 2014 4:29:18 AM (10 years ago)
Author:
Miguel Gonzalez <onzalezm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
beda555
Parents:
f44b076
Message:

Added the three models for hollow and massive rectangular parallelepipeds as defined in Z. Phys. Chem. 226 (2012) 837-854.
To do: Create model documentation.

Files:
6 added
4 edited

Legend:

Unmodified
Added
Removed

  • sansview/default_categories.json

    r27b7acc r201af9f  
    1 {"Shapes": [["PearlNecklaceModel", true], ["HollowCylinderModel", true], ["CoreMultiShellModel", true], ["FCCrystalModel", true], ["ParallelepipedModel", true], ["FuzzySphereModel", true], ["LamellarPCrystalModel", true], ["BCCrystalModel", true], ["LamellarFFHGModel", true], ["CSParallelepipedModel", true], ["OnionExpShellModel", true], ["MultiShellModel", true], ["RaspBerryModel", true], ["CoreShellEllipsoidModel", true], ["CappedCylinderModel", true], ["CylinderModel", true], ["SphericalSLDModel", true], ["SCCrystalModel", true], ["StackedDisksModel", true], ["LamellarPSHGModel", true], ["CoreShellBicelleModel", true], ["LamellarPSModel", true], ["EllipsoidModel", true], ["SphereModel", true], ["VesicleModel", true], ["TriaxialEllipsoidModel", true], ["LamellarModel", true], ["BarBellModel", true], ["CoreShellCylinderModel", true], ["BinaryHSModel", true], ["FlexibleCylinderModel", true], ["EllipticalCylinderModel", true], ["FlexCylEllipXModel", true], ["LinearPearlsModel", true], ["CoreShellModel", true]], "Structure Factor": [["SquareWellStructure", true], ["HayterMSAStructure", true], ["HardsphereStructure", true], ["StickyHSStructure", true]], "Shape-Independent": [["TwoPowerLawModel", true], ["GelFitModel", true], ["DABModel", true], ["PowerLawAbsModel", true], ["PorodModel", true], ["GuinierModel", true], ["TeubnerStreyModel", true], ["StarPolymer", true], ["BEPolyelectrolyte", true], ["UnifiedPowerRgModel", true], ["GuinierPorodModel", true], ["MassSurfaceFractal", true], ["Core2ndMomentModel", true], ["RPA10Model", true], ["TwoLorentzianModel", true], ["DebyeModel", true], ["MassFractalModel", true], ["Poly_GaussCoil", true], ["CorrLengthModel", true], ["PolymerExclVolume", true], ["FractalModel", true], ["PeakLorentzModel", true], ["BroadPeakModel", true], ["PeakGaussModel", true], ["FractalCoreShellModel", true], ["SurfaceFractalModel", true], ["GaussLorentzGelModel", true], ["LorentzModel", true]], "Uncategorized": [["ReflectivityModel", true], ["ReflectivityIIModel", true], ["PringlesModel", true], ["LineModel", true]]} 
     1{"Shapes": [["PearlNecklaceModel", true], ["HollowCylinderModel", true], ["CoreMultiShellModel", true], ["FCCrystalModel", true], ["ParallelepipedModel", true], ["FuzzySphereModel", true], ["LamellarPCrystalModel", true], ["BCCrystalModel", true], ["LamellarFFHGModel", true], ["CSParallelepipedModel", true], ["OnionExpShellModel", true], ["MultiShellModel", true], ["RaspBerryModel", true], ["CoreShellEllipsoidModel", true], ["CappedCylinderModel", true], ["CylinderModel", true], ["SphericalSLDModel", true], ["SCCrystalModel", true], ["StackedDisksModel", true], ["LamellarPSHGModel", true], ["CoreShellBicelleModel", true], ["LamellarPSModel", true], ["EllipsoidModel", true], ["SphereModel", true], ["VesicleModel", true], ["TriaxialEllipsoidModel", true], ["LamellarModel", true], ["BarBellModel", true], ["CoreShellCylinderModel", true], ["BinaryHSModel", true], ["FlexibleCylinderModel", true], ["EllipticalCylinderModel", true], ["FlexCylEllipXModel", true], ["LinearPearlsModel", true], ["CoreShellModel", true], ["RectangularPrismModel", true], ["RectangularHollowPrismModel", true],["RectangularHollowPrismInfThinWallsModel", true]], "Structure Factor": [["SquareWellStructure", true], ["HayterMSAStructure", true], ["HardsphereStructure", true], ["StickyHSStructure", true]], "Shape-Independent": [["TwoPowerLawModel", true], ["GelFitModel", true], ["DABModel", true], ["PowerLawAbsModel", true], ["PorodModel", true], ["GuinierModel", true], ["TeubnerStreyModel", true], ["StarPolymer", true], ["BEPolyelectrolyte", true], ["UnifiedPowerRgModel", true], ["GuinierPorodModel", true], ["MassSurfaceFractal", true], ["Core2ndMomentModel", true], ["RPA10Model", true], ["TwoLorentzianModel", true], ["DebyeModel", true], ["MassFractalModel", true], ["Poly_GaussCoil", true], ["CorrLengthModel", true], ["PolymerExclVolume", true], ["FractalModel", true], ["PeakLorentzModel", true], ["BroadPeakModel", true], ["PeakGaussModel", true], ["FractalCoreShellModel", true], ["SurfaceFractalModel", true], ["GaussLorentzGelModel", true], ["LorentzModel", true]], "Uncategorized": [["ReflectivityModel", true], ["ReflectivityIIModel", true], ["PringlesModel", true], ["LineModel", true]]} 

  • src/sans/models/c_extension/libigor/libCylinder.c

    r230f479 r201af9f  
    33103310} 
    33113311 
     3312 
     3313/* 
     3314 * Massive rectangular prism (a <= b <=c)  
     3315 * using eqns. (12)+(16) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012). 
     3316 * It is totally equivalent to the Parallelepiped model in this library, 
     3317 * but instead of using as input the a, b, and c sides of the prism, 
     3318 * it uses a and the ratios (b/a) and (c/a). 
     3319 * This allows to keep the shape when using the polydispersity. 
     3320 */ 
     3321 
     3322double 
     3323RectangularPrism(double dp[], double q) 
     3324{ 
     3325    int i, j; 
     3326    double scale, aa, bb, cc, delrho, bkg;               
     3327    int nordi=76;                               //order of integration 
     3328    int nordj=76; 
     3329    double Pi; 
     3330    double sumi, sumj; 
     3331    double v1a, v1b, v2a, v2b;          // upper and lower integration limits 
     3332    double answer; 
     3333    double theta, phi, arg, termA, termB, termC, AP, ahalf, bhalf, chalf; 
     3334    double vol, sldp, sld; 
     3335 
     3336    Pi = 4.0*atan(1.0); 
     3337 
     3338    scale = dp[0];       
     3339    aa = dp[1];          // parameter short_side 
     3340    bb = aa * dp[2];     // parameter b/a ratio 
     3341    cc = aa * dp[3];     // parameter c/a ratio 
     3342    sldp = dp[4];        // scattering length density of the object 
     3343    sld = dp[5];         // scattering length density of the solvent 
     3344    delrho = sldp - sld; 
     3345    bkg = dp[6]; 
     3346 
     3347    vol = aa*bb*cc; 
     3348    ahalf = 0.5 * aa; 
     3349    bhalf = 0.5 * bb; 
     3350    chalf = 0.5 * cc; 
     3351 
     3352   //Integration limits to use in Gaussian quadrature 
     3353    v1a = 0.; 
     3354    v1b = Pi/2.;  //theta integration limits 
     3355    v2a = 0.; 
     3356    v2b = Pi/2.;  //phi integration limits 
     3357 
     3358    sumi = 0.0; 
     3359    for(i=0;i<nordi;i++) { 
     3360 
     3361            theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;  //z values for gaussian quadrature 
     3362 
     3363            arg = q*chalf*cos(theta); 
     3364            if (fabs(arg) > 1.e-16) {termC = sin(arg) / arg;} else {termC = 1.0;}   
     3365 
     3366            sumj = 0.0; 
     3367            for(j=0;j<nordj;j++) { 
     3368 
     3369                phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0; //z values for gaussian quadrature 
     3370 
     3371                // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 
     3372 
     3373                arg = q*ahalf*sin(theta)*sin(phi);  
     3374                if (fabs(arg) > 1.e-16) {termA = sin(arg) / arg;} else {termA = 1.0;} 
     3375                
     3376                arg = q*bhalf*sin(theta)*cos(phi);  
     3377                if (fabs(arg) > 1.e-16) {termB = sin(arg) / arg;} else {termB = 1.0;}      
     3378                
     3379                AP = termA * termB * termC;   
     3380 
     3381                sumj += Gauss76Wt[j] * (AP*AP); 
     3382 
     3383            } 
     3384 
     3385            sumj = (v2b-v2a) * sumj / 2.0; 
     3386            sumi += Gauss76Wt[i] * sumj * sin(theta); 
     3387 
     3388    } 
     3389 
     3390    answer = (v1b-v1a)/2.0*sumi; 
     3391 
     3392    // Normalize by Pi (Eqn. 16).  
     3393    // The term (ABC)^2 does not appear because it was introduced before from the defs of termA, termB, termC 
     3394    // The factor 2 appears because the theta integral has been defined between 0 and pi/2, instead of 0 to pi 
     3395    answer = 2.0 * answer / Pi; //Form factor P(q) 
     3396 
     3397    // Multiply by contrast^2 
     3398    answer *= delrho*delrho; 
     3399 
     3400    //normalize by volume 
     3401    answer *= vol; 
     3402         
     3403    //convert to [cm-1] 
     3404    answer *= 1.0e8; 
     3405 
     3406    //Scale 
     3407    answer *= scale; 
     3408 
     3409    // add in the background 
     3410    answer += bkg; 
     3411 
     3412    return answer; 
     3413} 
     3414 
     3415 
     3416/* 
     3417 * Hollow rectangular prism (a <= b <=c) with infinitely thin walls 
     3418 * using eqn. (7-11) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012) 
     3419 */ 
     3420 
     3421double 
     3422RectangularHollowPrismInfThinWalls(double dp[], double q) 
     3423{ 
     3424    int i,j; 
     3425    double scale,aa,bb,cc,delrho,bkg; 
     3426    int nordi=76; 
     3427    int nordj=76; 
     3428    double Pi; 
     3429    double sumi, sumj; 
     3430    double v1a, v1b, v2a, v2b;          //upper and lower integration limits 
     3431    double answer; 
     3432    double theta, phi, termAL_theta, termAT_theta, AL, AT, ahalf, bhalf, chalf, vol; 
     3433    double sldp, sld ; 
     3434 
     3435    Pi = 4.0*atan(1.0); 
     3436 
     3437    scale = dp[0];                //make local copies in case memory moves 
     3438    aa = dp[1]; 
     3439    bb = aa * dp[2]; 
     3440    cc = aa * dp[3]; 
     3441    sldp = dp[4];        // scattering length density of the object 
     3442    sld = dp[5];         // scattering length density of the solvent 
     3443    delrho = sldp - sld; 
     3444    bkg = dp[6]; 
     3445 
     3446    ahalf = 0.5 * aa; 
     3447    bhalf = 0.5 * bb; 
     3448    chalf = 0.5 * cc; 
     3449 
     3450    //Integration limits to use in Gaussian quadrature 
     3451    v1a = 0.; 
     3452    v1b = Pi/2;   //theta integration limits 
     3453    v2a = 0.; 
     3454    v2b = Pi/2.;  //phi integration limits 
     3455 
     3456    sumi = 0.0; 
     3457    for(i=0;i<nordi;i++) { 
     3458 
     3459        theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;      //z values for gaussian quadrature 
     3460 
     3461        // To check potential problems if denominator goes to zero here !!! 
     3462        termAL_theta = 8.0*cos(q*chalf*cos(theta)) / (q*q*sin(theta)*sin(theta)); 
     3463        termAT_theta = 8.0*sin(q*chalf*cos(theta)) / (q*q*sin(theta)*cos(theta)); 
     3464 
     3465        sumj = 0.0; 
     3466        for(j=0;j<nordj;j++) { 
     3467 
     3468            phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0;    //z values for gaussian quadrature 
     3469 
     3470            // Amplitude AL from eqn. (7c) 
     3471            AL = termAL_theta * sin(q*ahalf*sin(theta)*sin(phi)) * sin(q*bhalf*sin(theta)*cos(phi)) / (sin(phi)*cos(phi)); 
     3472 
     3473            // Amplitude AT from eqn. (9) 
     3474            AT = termAT_theta * (  (cos(q*ahalf*sin(theta)*sin(phi))*sin(q*bhalf*sin(theta)*cos(phi))/cos(phi)) + (cos(q*bhalf*sin(theta)*cos(phi))*sin(q*ahalf*sin(theta)*sin(phi))/sin(phi)) ); 
     3475 
     3476            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT); 
     3477 
     3478        } 
     3479 
     3480        sumj = (v2b-v2a) * sumj / 2.0; 
     3481        sumi += Gauss76Wt[i] * sumj * sin(theta); 
     3482 
     3483    } 
     3484 
     3485    answer = (v1b-v1a)/2.0*sumi; 
     3486 
     3487    // Normalize as in Eqn. (11). The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
     3488    vol = 2.0*aa*bb+2.0*aa*cc+2.0*bb*cc; 
     3489    answer = 2.0 * answer / Pi / vol / vol; 
     3490 
     3491    // Multiply by contrast^2 
     3492    answer *= delrho*delrho; 
     3493 
     3494    //normalize by volume 
     3495    answer *= vol; 
     3496 
     3497    //convert to [cm-1] 
     3498    answer *= 1.0e8; 
     3499 
     3500    //Scale 
     3501    answer *= scale; 
     3502 
     3503    // add in the background 
     3504    answer += bkg; 
     3505 
     3506    return answer; 
     3507} 
     3508 
     3509 
     3510/* 
     3511 * Hollow rectangular prism (a <= b <=c) 
     3512 * using eqn. (13-15) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012) 
     3513 */ 
     3514  
     3515double 
     3516RectangularHollowPrism(double dp[], double q) 
     3517{ 
     3518    int i, j; 
     3519    double scale, aa, bb, cc, thickness, delrho, bkg; 
     3520    int nordi=76; 
     3521    int nordj=76; 
     3522    double Pi; 
     3523    double sumi, sumj; 
     3524    double v1a, v1b, v2a, v2b;          //upper and lower integration limits 
     3525    double answer; 
     3526    double theta, phi, arg, termA1, termB1, termC1, termA2, termB2, termC2, AP1, AP2, AP, ahalf, bhalf, chalf, vol; 
     3527    double sldp, sld ; 
     3528 
     3529    Pi = 4.0*atan(1.0); 
     3530 
     3531    scale = dp[0]; 
     3532    aa = dp[1]; 
     3533    bb = aa * dp[2]; 
     3534    cc = aa * dp[3]; 
     3535    thickness = dp[4]; 
     3536    sldp = dp[5];        // scattering length density of the object 
     3537    sld = dp[6];         // scattering length density of the solvent 
     3538    delrho = sldp - sld; 
     3539    bkg = dp[7]; 
     3540 
     3541    ahalf = 0.5 * aa; 
     3542    bhalf = 0.5 * bb; 
     3543    chalf = 0.5 * cc; 
     3544 
     3545    //Integration limits to use in Gaussian quadrature 
     3546    v1a = 0.; 
     3547    v1b = Pi/2;    //theta integration limits 
     3548    v2a = 0.; 
     3549    v2b = Pi/2.;   //phi integration limits 
     3550 
     3551    sumi = 0.0; 
     3552    for(i=0;i<nordi;i++) { 
     3553 
     3554        theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;      //z values for gaussian quadrature 
     3555 
     3556            arg = q*chalf*cos(theta); 
     3557            if (fabs(arg) > 1.e-16) {termC1 = sin(arg) / arg;} else {termC1 = 1.0;} 
     3558            arg = q*(chalf-thickness)*cos(theta); 
     3559            if (fabs(arg) > 1.e-16) {termC2 = sin(arg) / arg;} else {termC2 = 1.0;} 
     3560 
     3561            sumj = 0.0; 
     3562            for(j=0;j<nordj;j++) { 
     3563 
     3564            phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0;    //z values for gaussian quadrature 
     3565 
     3566            // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0 
     3567 
     3568                arg = q*ahalf*sin(theta)*sin(phi); 
     3569                if (fabs(arg) > 1.e-16) {termA1 = sin(arg) / arg;} else {termA1 = 1.0;} 
     3570                arg = q*(ahalf-thickness)*sin(theta)*sin(phi); 
     3571                if (fabs(arg) > 1.e-16) {termA2 = sin(arg) / arg;} else {termA2 = 1.0;} 
     3572 
     3573                arg = q*bhalf*sin(theta)*cos(phi); 
     3574                if (fabs(arg) > 1.e-16) {termB1 = sin(arg) / arg;} else {termB1 = 1.0;} 
     3575                arg = q*(bhalf-thickness)*sin(theta)*cos(phi); 
     3576                if (fabs(arg) > 1.e-16) {termB2 = sin(arg) / arg;} else {termB2 = 1.0;} 
     3577 
     3578            AP1 = (aa*bb*cc) * termA1 * termB1 * termC1; 
     3579            AP2 = 8 * (ahalf-thickness) * (bhalf-thickness) * (chalf-thickness) * termA2 * termB2 * termC2; 
     3580            AP = AP1 - AP2; 
     3581 
     3582            sumj += Gauss76Wt[j] * (AP*AP); 
     3583 
     3584            } 
     3585 
     3586            sumj = (v2b-v2a) * sumj / 2.0; 
     3587            sumi += Gauss76Wt[i] * sumj * sin(theta); 
     3588 
     3589    } 
     3590 
     3591    answer = (v1b-v1a)/2.0*sumi; 
     3592 
     3593    // Normalize as in Eqn. (15). 
     3594    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
     3595    vol = (aa*bb*cc-((aa-2*thickness)*(bb-2*thickness)*(cc-2*thickness))); 
     3596    answer = 2.0 * answer / Pi / vol / vol; 
     3597 
     3598    // Multiply by contrast^2 
     3599    answer *= delrho*delrho; 
     3600 
     3601    //normalize by volume 
     3602    answer *= vol; 
     3603 
     3604    //convert to [cm-1] 
     3605    answer *= 1.0e8; 
     3606 
     3607    //Scale 
     3608    answer *= scale; 
     3609 
     3610    // add in the background 
     3611    answer += bkg; 
     3612 
     3613    return answer; 
     3614} 
     3615 
     3616 

  • src/sans/models/c_extension/libigor/libCylinder.h

    r230f479 r201af9f  
    3838double PolyCoreBicelle(double dp[], double q); 
    3939double CSParallelepiped(double dp[], double q); 
     40double RectangularPrism(double dp[], double q); 
     41double RectangularHollowPrismInfThinWalls(double dp[], double q); 
     42double RectangularHollowPrism(double dp[], double q); 
    4043 
    4144/* internal functions */ 

  • src/sans/perspectives/fitting/models.py

    r27b7acc r201af9f  
    677677        self.model_name_list.append(PringlesModel.__name__) 
    678678 
     679        from sans.models.RectangularPrismModel import RectangularPrismModel 
     680        self.model_dictionary[RectangularPrismModel.__name__] = RectangularPrismModel 
     681        self.shape_list.append(RectangularPrismModel) 
     682        self.multiplication_factor.append(RectangularPrismModel) 
     683        self.model_name_list.append(RectangularPrismModel.__name__) 
     684 
     685        from sans.models.RectangularHollowPrismInfThinWallsModel import RectangularHollowPrismInfThinWallsModel 
     686        self.model_dictionary[RectangularHollowPrismInfThinWallsModel.__name__] = RectangularHollowPrismInfThinWallsModel 
     687        self.shape_list.append(RectangularHollowPrismInfThinWallsModel) 
     688        self.multiplication_factor.append(RectangularHollowPrismInfThinWallsModel) 
     689        self.model_name_list.append(RectangularHollowPrismInfThinWallsModel.__name__) 
     690 
     691        from sans.models.RectangularHollowPrismModel import RectangularHollowPrismModel 
     692        self.model_dictionary[RectangularHollowPrismModel.__name__] = RectangularHollowPrismModel 
     693        self.shape_list.append(RectangularHollowPrismModel) 
     694        self.multiplication_factor.append(RectangularHollowPrismModel) 
     695        self.model_name_list.append(RectangularHollowPrismModel.__name__) 
     696 
    679697        #from sans.models.FractalO_Z import FractalO_Z 
    680698        #self.model_dictionary[FractalO_Z.__name__] = FractalO_Z 
Note: See TracChangeset for help on using the changeset viewer.