Changeset 633 for sans/Dev


Ignore:
Timestamp:
Mar 3, 2010 9:05:47 AM (13 years ago)
Author:
srkline
Message:

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.

Location:
sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models
Files:
28 edited

Legend:

Unmodified
Added
Removed

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/CoreShellCylinder_v40.ipf

    r570 r633  
    205205   //Local variables  
    206206        Variable dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval 
     207        Variable si1,si2,be1,be2 
    207208         
    208209        dr1 = rhoc-rhos 
    … …  
    216217        sinarg2 = qq*(length+thick)*cos(dum) 
    217218         
    218         t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1 
    219         t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2 
     219        if(besarg1 == 0.0) 
     220                be1 = 0.5 
     221        else 
     222                be1 = bessJ(1,besarg1)/besarg1 
     223        endif 
     224         
     225        if(besarg2 == 0.0) 
     226                be2 = 0.5 
     227        else 
     228                be2 = bessJ(1,besarg2)/besarg2 
     229        endif 
     230         
     231        if(sinarg1 == 0.0) 
     232                si1 = 1.0 
     233        else 
     234                si1 = sin(sinarg1)/sinarg1 
     235        endif 
     236         
     237        if(sinarg2 == 0.0) 
     238                si2 = 1.0 
     239        else 
     240                si2 = sin(sinarg2)/sinarg2 
     241        endif 
     242 
     243        t1 = 2.0*vol1*dr1*si1*be1 
     244        t2 = 2.0*vol2*dr2*si2*be2 
    220245         
    221246        retval = ((t1+t2)^2)*sin(dum) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/CoreShell_v40.ipf

    r570 r633  
    120120        qr=x*rcore 
    121121        contr = rhocore-rhoshel 
    122         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     122        if(qr==0) 
     123                bes = 1 
     124        else 
     125                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     126        endif 
    123127        vol = 4*pi/3*rcore^3 
    124128        f = vol*bes*contr 
    125129        //now the shell 
    126130        qr=x*(rcore+thick) 
     131        if(qr==0) 
     132                bes = 1 
     133        else 
     134                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     135        endif 
    127136        contr = rhoshel-rhosolv 
    128         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
    129137        vol = 4*pi/3*(rcore+thick)^3 
    130138        f += vol*bes*contr 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/Cylinder_v40.ipf

    r583 r633  
    197197 
    198198   //Local variables  
    199         Variable besarg,bj,retval,d1,t1,b1,t2,b2 
     199        Variable besarg,bj,retval,d1,t1,b1,t2,b2,be,si,siarg 
    200200    besarg = qq*rr*sin(theta) 
    201  
     201        siarg = qq * h * cos(theta) 
     202         
    202203    bj =bessJ(1,besarg) 
    203204 
    204205//* Computing 2nd power */ 
    205     d1 = sin(qq * h * cos(theta)) 
     206    d1 = sin(siarg) 
    206207    t1 = d1 * d1 
    207208//* Computing 2nd power */ 
    … …  
    209210    t2 = d1 * d1 * 4.0 * sin(theta) 
    210211//* Computing 2nd power */ 
    211     d1 = qq * h * cos(theta) 
     212    d1 = siarg 
    212213    b1 = d1 * d1 
    213214//* Computing 2nd power */ 
    214215    d1 = qq * rr * sin(theta) 
    215216    b2 = d1 * d1 
    216     retval = t1 * t2 / b1 / b2 
     217     
     218    if(besarg == 0) 
     219        be = sin(theta) 
     220    else 
     221        be = t2/b2 
     222    endif 
     223     
     224    if(siarg == 0) 
     225        si = 1 
     226    else 
     227        si = t1/b1 
     228    endif 
     229     
     230    retval = be*si 
    217231 
    218232    return retval 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/HardSphereStruct_v40.ipf

    r570 r633  
    6868        //[0] radius 
    6969        //[1] volume fraction 
    70         //Variable timer=StartMSTimer 
    7170         
    7271        Variable r,phi,struc 
    … …  
    105104        STRUC = VSTRUC 
    106105//   10 CONTINUE 
    107         //Variable elapse=StopMSTimer(timer) 
    108       //Print "HS struct eval time (s) = ",elapse 
     106 
    109107      RETURN Struc 
    110108End 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/HollowCylinders_v40.ipf

    r570 r633  
    193193 
    194194   //Local variables  
    195         Variable gamma,besarg1,besarg2,lam1,lam2,t2,retval,psi,sinarg 
    196          
    197         gamma = r2/r1 
     195        Variable gamm,besarg1,besarg2,lam1,lam2,t2,retval,psi,sinarg 
     196         
     197        gamm = r2/r1 
    198198        besarg1 = qq*r1*sqrt(1-theta^2) 
    199199        besarg2 = qq*r2*sqrt(1-theta^2) 
    200         lam1 = 2*bessJ(1,besarg1)/besarg1 
    201         lam2 = 2*bessJ(1,besarg2)/besarg2 
    202         psi = 1/(1-gamma^2)*(lam1 -  gamma^2*lam2)              //SRK 10/19/00 
     200         
     201        if(besarg1 == 0) 
     202                lam1 = 1 
     203        else 
     204                lam1 = 2*bessJ(1,besarg1)/besarg1 
     205        endif 
     206        if(besarg2 == 0) 
     207                lam2 = 1 
     208        else 
     209                lam2 = 2*bessJ(1,besarg2)/besarg2 
     210        endif 
     211         
     212        psi = 1/(1-gamm^2)*(lam1 -  gamm^2*lam2)                //SRK 10/19/00 
    203213         
    204214        sinarg = qq*h*theta/2 
    205         t2 = sin(sinarg)/sinarg 
    206          
     215        if(sinarg == 0) 
     216                t2 = 1 
     217        else 
     218                t2 = sin(sinarg)/sinarg 
     219        endif 
     220                 
    207221        retval = psi*psi*t2*t2 
    208222         

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Beaucage_v40.ipf

    r570 r633  
    285285        ans = G1*exp(-x*x*Rg1*Rg1/3) 
    286286        ans += B1*(erf1^3/x)^Pow1 
     287         
     288        if(x == 0) 
     289                ans = G1 
     290        endif 
     291         
    287292        ans *= scale 
    288293        ans += bkg 
    289294         
    290         Return (ans) 
     295        return(ans) 
    291296End 
    292297 
    … …  
    332337        ans += G2*exp(-x*x*Rg2*Rg2/3) 
    333338        ans += B2*(erf2^3/x)^Pow2 
     339         
     340        if(x == 0) 
     341                ans = G1 + G2 
     342        endif 
     343         
    334344        ans *= scale 
    335345        ans += bkg 
    … …  
    384394        ans += G2*exp(-x*x*Rg2*Rg2/3) + B2*exp(-x*x*Rg3*Rg3/3)*(erf2^3/x)^Pow2 
    385395        ans += G3*exp(-x*x*Rg3*Rg3/3) + B3*(erf3^3/x)^Pow3 
     396         
     397        if(x == 0) 
     398                ans = G1 + G2 + G3 
     399        endif 
     400         
    386401        ans *= scale 
    387402        ans += bkg 
    … …  
    442457        ans += G3*exp(-x*x*Rg3*Rg3/3) + B3*exp(-x*x*Rg4*Rg4/3)*(erf3^3/x)^Pow3 
    443458        ans += G4*exp(-x*x*Rg4*Rg4/3) + B4*(erf4^3/x)^Pow4 
     459         
     460        if(x == 0) 
     461                ans = G1 + G2 + G3 + G4 
     462        endif 
     463         
    444464        ans *= scale 
    445465        ans += bkg 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Cylinder_PolyLength_v40.ipf

    r570 r633  
    203203        //calculate the orientationally averaged P(q) for the input rad 
    204204        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170. 
    205         Make/O/n=6 kernpar 
     205        Make/O/D/n=6 kernpar 
    206206        Wave kp = kernpar 
    207207        kp[0] = 1               //scale fixed at 1 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/FlexCyl_PolyLen_v40.ipf

    r570 r633  
    197197        //ww[6] = bkg [cm-1] 
    198198        Variable Pq,vcyl,dl 
    199         Make/O/n=7 fle_ker 
     199        Make/O/D/n=7 fle_ker 
    200200        Wave kp = fle_ker 
    201201        kp[0] = 1               //scale fixed at 1 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/FlexCyl_PolyRadius_v40.ipf

    r570 r633  
    201201        //calculate the orientationally averaged P(q) for the input rad 
    202202        //this is correct - see K&C (1983) or Lin &Tsao JACryst (1996)29 170. 
    203         Make/O/n=7 kernpar 
     203        Make/O/D/n=7 kernpar 
    204204        Wave kp = kernpar 
    205205        kp[0] = 1               //scale fixed at 1 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/FlexibleCylinder_v40.ipf

    r570 r633  
    121121        Variable flex,crossSect 
    122122        flex = Sk_WR(x,L,B) 
    123       
    124         crossSect = (2*bessJ(1,qr)/qr)^2 
     123     
     124    if(qr == 0) 
     125            crossSect = 1 
     126    else 
     127        crossSect = (2*bessJ(1,qr)/qr)^2 
     128    endif 
    125129         
    126130        //normalize form factor by multiplying by cylinder volume * cont^2 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Fractal_v40.ipf

    r570 r633  
    121121        variable x 
    122122         
    123         variable r0,Df,corr,phi,sldp,sldm,bkg 
     123        variable r0,Df,corr,phi,sldp,sldm,bkg,bes 
    124124        variable pq,sq,ans 
    125125        phi=w[0]   // volume fraction of building block spheres... 
    … …  
    132132         
    133133        //calculate P(q) for the spherical subunits, units cm-1 sr-1 
    134         pq = 1.0e8*phi*4/3*pi*r0^3*(sldp-sldm)^2*(3*(sin(x*r0) - x*r0*cos(x*r0))/(x*r0)^3)^2 
     134        if(x*r0 == 0) 
     135                bes = 1 
     136        else 
     137                bes = (3*(sin(x*r0) - x*r0*cos(x*r0))/(x*r0)^3)^2 
     138        endif 
     139         
     140        pq = 1.0e8*phi*4/3*pi*r0^3*(sldp-sldm)^2*bes 
    135141 
    136142        //calculate S(q) 
    … …  
    138144        sq /= (x*r0)^Df * (1 + 1/(x*corr)^2)^((Df-1)/2) 
    139145        sq += 1 
     146 
    140147        //combine and return 
    141148        ans = pq*sq + bkg 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/MultiShell_v40.ipf

    r570 r633  
    141141         
    142142        Variable val=0 
    143          
    144         val = 3*(sin(qr) - qr*cos(qr))/qr^3 
    145          
     143        if(qr == 0) 
     144                val = 1 
     145        else     
     146                val = 3*(sin(qr) - qr*cos(qr))/qr^3 
     147        endif 
     148                 
    146149        return(val) 
    147150End 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/PolyCoreShellCylinder_v40.ipf

    r570 r633  
    285285 
    286286        Variable dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval         //Local variables  
    287  
     287        Variable si1,si2,be1,be2 
     288         
    288289        dr1 = rhoc-rhos 
    289290        dr2 = rhos-rhosolv 
    … …  
    296297        sinarg2 = qq*(length+facthick)*cos(dum) 
    297298         
    298         t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1 
    299         t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2 
     299        if(besarg1 == 0) 
     300                be1 = 0.5 
     301        else 
     302                be1 = bessJ(1,besarg1)/besarg1 
     303        endif 
     304        if(besarg2 == 0) 
     305                be2 = 0.5 
     306        else 
     307                be2 = bessJ(1,besarg2)/besarg2 
     308        endif 
     309        if(sinarg1 == 0) 
     310                si1 = 1 
     311        else 
     312                si1 = sin(sinarg1)/sinarg1 
     313        endif 
     314        if(sinarg2 == 0) 
     315                si2 = 1 
     316        else 
     317                si2 = sin(sinarg2)/sinarg2 
     318        endif 
     319         
     320         
     321        t1 = 2*vol1*dr1*si1*be1 
     322        t2 = 2*vol2*dr2*si2*be2 
    300323         
    301324        retval = ((t1+t2)^2)*sin(dum) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/TriaxialEllipsoid_v40.ipf

    r570 r633  
    186186        arg = qq*sqrt(arg) 
    187187         
    188         val = 9*((sin(arg) - arg*cos(arg))/arg^3 )^2 
    189          
     188        if(arg == 0) 
     189                val = 1 
     190        else 
     191                val = 9*((sin(arg) - arg*cos(arg))/arg^3 )^2 
     192        endif 
     193 
    190194        return(val) 
    191195end 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/Vesicle_UL_v40.ipf

    r570 r633  
    131131        qr=x*rcore 
    132132        contr = rhocore-rhoshel 
    133         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     133        if(qr == 0) 
     134                bes = 1 
     135        else 
     136                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     137        endif 
    134138        vol = 4*pi/3*rcore^3 
    135139        f = vol*bes*contr 
    … …  
    137141        qr=x*(rcore+thick) 
    138142        contr = rhoshel-rhosolv 
    139         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     143         
     144        if(qr == 0) 
     145                bes = 1 
     146        else 
     147                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     148        endif 
     149         
    140150        vol = 4*pi/3*(rcore+thick)^3 
    141151        f += vol*bes*contr 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Barbell_v40.ipf

    r570 r633  
    149149        Variable retVal 
    150150        Variable scale,contr,bkg,inten,sldc,slds 
    151         Variable len,rad,hDist,endRad 
     151        Variable len,rad,hDist,endRad,be 
    152152        scale = w[0] 
    153153        rad = w[1] 
    … …  
    169169        arg2 = x*rad*sin(dum) 
    170170         
    171         retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2 
     171        if(arg2 == 0) 
     172                be = 0.5 
     173        else 
     174                be = Besselj(1, arg2)/arg2 
     175        endif 
     176        retVal += pi*rad*rad*len*sinc(arg1)*2*be 
    172177         
    173178        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta) 
    … …  
    209214        Variable val,arg1,arg2 
    210215        Variable scale,contr,bkg,inten,sldc,slds 
    211         Variable len,rad,hDist,endRad 
     216        Variable len,rad,hDist,endRad,be 
    212217        scale = w[0] 
    213218        rad = w[1] 
    … …  
    223228        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt) 
    224229         
    225         val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2 
     230        if(arg2 == 0) 
     231                be = 0.5 
     232        else 
     233                be = Besselj(1, arg2)/arg2 
     234        endif 
     235        val = cos(arg1)*(1-tt*tt)*be 
    226236         
    227237        return(val) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/CappedCylinder_v40.ipf

    r570 r633  
    167167        retval = IntegrateFn76(CapCyl_Inner,-hDist/endRad,1,w,x) 
    168168         
    169         Variable arg1,arg2 
     169        Variable arg1,arg2,be 
    170170        arg1 = x*len/2*cos(dum) 
    171171        arg2 = x*rad*sin(dum) 
    172172         
    173         retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2 
     173        if(arg2 == 0) 
     174                be = 0.5 
     175        else 
     176                be = Besselj(1, arg2)/arg2 
     177        endif 
     178        retVal += pi*rad*rad*len*sinc(arg1)*2*be 
    174179         
    175180        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta) 
    … …  
    211216        Variable val,arg1,arg2 
    212217        Variable scale,contr,bkg,inten,sldc,slds 
    213         Variable len,rad,hDist,endRad 
     218        Variable len,rad,hDist,endRad,be 
    214219        scale = w[0] 
    215220        rad = w[1] 
    … …  
    225230        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt) 
    226231         
    227         val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2 
     232        if(arg2 == 0) 
     233                be = 0.5 
     234        else 
     235                be = Besselj(1, arg2)/arg2 
     236        endif 
     237         
     238        val = cos(arg1)*(1-tt*tt)*be 
    228239         
    229240        return(val) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/ConvexLens_v40.ipf

    r570 r633  
    179179        retval = IntegrateFn76(ConvLens_Inner,-hDist/endRad,1,w,x) 
    180180         
    181         Variable arg1,arg2 
     181        Variable arg1,arg2,be 
    182182        arg1 = x*len/2*cos(dum) 
    183183        arg2 = x*rad*sin(dum) 
    184184         
    185         retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2 
     185        if(arg2 == 0) 
     186                be = 0.5 
     187        else 
     188                be = Besselj(1, arg2)/arg2 
     189        endif 
     190         
     191        retVal += pi*rad*rad*len*sinc(arg1)*2*be 
    186192         
    187193        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta) 
    … …  
    223229        Variable val,arg1,arg2 
    224230        Variable scale,contr,bkg,inten,sldc,slds 
    225         Variable len,rad,hDist,endRad 
     231        Variable len,rad,hDist,endRad,be 
    226232        scale = w[0] 
    227233        rad = w[1] 
    … …  
    237243        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt) 
    238244         
    239         val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2 
     245        if(arg2 == 0) 
     246                be = 0.5 
     247        else 
     248                be = Besselj(1, arg2)/arg2 
     249        endif 
     250         
     251        val = cos(arg1)*(1-tt*tt)*be 
    240252         
    241253        return(val) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Core_and_NShells_v40.ipf

    r570 r633  
    363363        qr=x*rcore 
    364364        contr = rhocore-rhoshel 
    365         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     365         
     366        if(qr == 0) 
     367                bes = 1 
     368        else 
     369                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     370        endif 
     371         
    366372        vol = 4*pi/3*rcore^3 
    367373        f = vol*bes*contr 
    … …  
    369375        qr=x*(rcore+thick) 
    370376        contr = rhoshel-rhosolv 
    371         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     377         
     378        if(qr == 0) 
     379                bes = 1 
     380        else 
     381                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     382        endif 
     383         
    372384        vol = 4*pi/3*(rcore+thick)^3 
    373385        f += vol*bes*contr 
    … …  
    421433        qr=x*rcore 
    422434        contr = rhocore-rhoshel1 
    423         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     435        if(qr == 0) 
     436                bes = 1 
     437        else 
     438                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     439        endif 
    424440        vol = 4*pi/3*rcore^3 
    425441        f = vol*bes*contr 
    … …  
    427443        qr=x*(rcore+thick1) 
    428444        contr = rhoshel1-rhoshel2 
    429         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     445        if(qr == 0) 
     446                bes = 1 
     447        else 
     448                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     449        endif 
    430450        vol = 4*pi/3*(rcore+thick1)^3 
    431451        f += vol*bes*contr 
    … …  
    433453        qr=x*(rcore+thick1+thick2) 
    434454        contr = rhoshel2-rhosolv 
    435         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     455        if(qr == 0) 
     456                bes = 1 
     457        else 
     458                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     459        endif 
    436460        vol = 4*pi/3*(rcore+thick1+thick2)^3 
    437461        f += vol*bes*contr 
    … …  
    488512        qr=x*rcore 
    489513        contr = rhocore-rhoshel1 
    490         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     514        if(qr == 0) 
     515                bes = 1 
     516        else 
     517                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     518        endif 
    491519        vol = 4*pi/3*rcore^3 
    492520        f = vol*bes*contr 
    … …  
    494522        qr=x*(rcore+thick1) 
    495523        contr = rhoshel1-rhoshel2 
    496         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     524        if(qr == 0) 
     525                bes = 1 
     526        else 
     527                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     528        endif 
    497529        vol = 4*pi/3*(rcore+thick1)^3 
    498530        f += vol*bes*contr 
    … …  
    500532        qr=x*(rcore+thick1+thick2) 
    501533        contr = rhoshel2-rhoshel3 
    502         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     534        if(qr == 0) 
     535                bes = 1 
     536        else 
     537                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     538        endif 
    503539        vol = 4*pi/3*(rcore+thick1+thick2)^3 
    504540        f += vol*bes*contr 
    … …  
    506542        qr=x*(rcore+thick1+thick2+thick3) 
    507543        contr = rhoshel3-rhosolv 
    508         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     544        if(qr == 0) 
     545                bes = 1 
     546        else 
     547                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     548        endif 
    509549        vol = 4*pi/3*(rcore+thick1+thick2+thick3)^3 
    510550        f += vol*bes*contr 
    … …  
    564604        qr=x*rcore 
    565605        contr = rhocore-rhoshel1 
    566         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     606        if(qr == 0) 
     607                bes = 1 
     608        else 
     609                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     610        endif 
    567611        vol = 4*pi/3*rcore^3 
    568612        f = vol*bes*contr 
    … …  
    570614        qr=x*(rcore+thick1) 
    571615        contr = rhoshel1-rhoshel2 
    572         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     616        if(qr == 0) 
     617                bes = 1 
     618        else 
     619                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     620        endif 
    573621        vol = 4*pi/3*(rcore+thick1)^3 
    574622        f += vol*bes*contr 
    … …  
    576624        qr=x*(rcore+thick1+thick2) 
    577625        contr = rhoshel2-rhoshel3 
    578         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     626        if(qr == 0) 
     627                bes = 1 
     628        else 
     629                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     630        endif 
    579631        vol = 4*pi/3*(rcore+thick1+thick2)^3 
    580632        f += vol*bes*contr 
    … …  
    582634        qr=x*(rcore+thick1+thick2+thick3) 
    583635        contr = rhoshel3-rhoshel4 
    584         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     636        if(qr == 0) 
     637                bes = 1 
     638        else 
     639                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     640        endif 
    585641        vol = 4*pi/3*(rcore+thick1+thick2+thick3)^3 
    586642        f += vol*bes*contr 
    … …  
    588644        qr=x*(rcore+thick1+thick2+thick3+thick4) 
    589645        contr = rhoshel4-rhosolv 
    590         bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     646        if(qr == 0) 
     647                bes = 1 
     648        else 
     649                bes = 3*(sin(qr)-qr*cos(qr))/qr^3 
     650        endif 
    591651        vol = 4*pi/3*(rcore+thick1+thick2+thick3+thick4)^3 
    592652        f += vol*bes*contr 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Dumbbell_v40.ipf

    r570 r633  
    162162        Variable retVal 
    163163        Variable scale,contr,bkg,inten,sldc,slds 
    164         Variable len,rad,hDist,endRad 
     164        Variable len,rad,hDist,endRad,be 
    165165        scale = w[0] 
    166166        rad = w[1] 
    … …  
    182182        arg2 = x*rad*sin(dum) 
    183183         
    184         retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2 
     184        if(arg2 == 0) 
     185                be = 0.5 
     186        else 
     187                be = Besselj(1, arg2)/arg2 
     188        endif 
     189         
     190        retVal += pi*rad*rad*len*sinc(arg1)*2*be 
    185191         
    186192        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta) 
    … …  
    223229        Variable val,arg1,arg2 
    224230        Variable scale,contr,bkg,inten,sldc,slds 
    225         Variable len,rad,hDist,endRad 
     231        Variable len,rad,hDist,endRad,be 
    226232        scale = w[0] 
    227233        rad = w[1] 
    … …  
    237243        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt) 
    238244         
    239         val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2 
     245        if(arg2 == 0) 
     246                be = 0.5 
     247        else 
     248                be = Besselj(1, arg2)/arg2 
     249        endif 
     250        val = cos(arg1)*(1-tt*tt)*be 
    240251         
    241252        return(val) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/PolyGaussShell_v40.ipf

    r570 r633  
    155155         
    156156        retval += bkg 
    157          
     157                 
    158158        return(retval) 
    159159         

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Spherocylinder_v40.ipf

    r570 r633  
    124124        bkg = w[5] 
    125125         
     126        endRad = rad 
     127         
    126128        Make/O/D/N=7 SphCyl_tmp 
    127129        SphCyl_tmp[0] = w[0] 
    … …  
    158160        Variable retVal 
    159161        Variable scale,contr,bkg,inten,sldc,slds 
    160         Variable len,rad,hDist,endRad 
     162        Variable len,rad,hDist,endRad,be 
    161163        scale = w[0] 
    162164        rad = w[1] 
    … …  
    178180        arg2 = x*rad*sin(dum) 
    179181         
    180         retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2 
     182        if(arg2 == 0) 
     183                be = 0.5 
     184        else 
     185                be = Besselj(1, arg2)/arg2 
     186        endif 
     187         
     188        retVal += pi*rad*rad*len*sinc(arg1)*2*be 
    181189         
    182190        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta) 
    … …  
    218226        Variable val,arg1,arg2 
    219227        Variable scale,contr,bkg,inten,sldc,slds 
    220         Variable len,rad,hDist,endRad 
     228        Variable len,rad,hDist,endRad,be 
    221229        scale = w[0] 
    222230        rad = w[1] 
    … …  
    232240        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt) 
    233241         
    234         val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2 
     242        if(arg2 == 0) 
     243                be = 0.5 
     244        else 
     245                be = Besselj(1, arg2)/arg2 
     246        endif 
     247         
     248        val = cos(arg1)*(1-tt*tt)*be 
    235249         
    236250        return(val) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2009/PolyCoreBicelle_v40.ipf

    r570 r633  
    269269 
    270270        Variable dr1,dr2, dr3, besarg1,besarg2, vol1,vol2, vol3,sinarg1,sinarg2,t1,t2,t3, retval                //Local variables  
     271        Variable si1,si2,be1,be2 
    271272 
    272273        dr1 = rhoc-rhoh 
    … …  
    281282        sinarg2 = qq*(length+facthick)*cos(dum) 
    282283         
    283         t1 = 2*vol1*dr1*sin(sinarg1)/sinarg1*bessJ(1,besarg1)/besarg1 
    284         t2 = 2*vol2*dr2*sin(sinarg2)/sinarg2*bessJ(1,besarg2)/besarg2 
    285         t3 = 2*vol3*dr3*sin(sinarg2)/sinarg2*bessJ(1,besarg1)/besarg1 
     284        if(besarg1 == 0) 
     285                be1 = 0.5 
     286        else 
     287                be1 = bessJ(1,besarg1)/besarg1 
     288        endif 
     289        if(besarg2 == 0) 
     290                be2 = 0.5 
     291        else 
     292                be2 = bessJ(1,besarg2)/besarg2 
     293        endif 
     294        if(sinarg1 == 0) 
     295                si1 = 1 
     296        else 
     297                si1 = sin(sinarg1)/sinarg1 
     298        endif 
     299        if(sinarg2 == 0) 
     300                si2 = 1 
     301        else 
     302                si2 = sin(sinarg2)/sinarg2 
     303        endif 
     304         
     305        t1 = 2*vol1*dr1*si1*be1 
     306        t2 = 2*vol2*dr2*si2*be2 
     307        t3 = 2*vol3*dr3*si2*be1 
    286308         
    287309        retval = ((t1+t2+t3)^2)*sin(dum) 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2010/PolymerExcludVol_v40.ipf

    r619 r633  
    121121        Ps=(1/(nu*Xx^o2nu))*(gammaInc(o2nu,Xx,0)-(1/Xx^o2nu)*gammaInc(onu,Xx,0)) 
    122122        Debye=(2/Xx^2)*(exp(-Xx)-1+Xx) 
     123         
     124        if(qval == 0) 
     125                Ps = 1 
     126        endif 
     127         
    123128        inten = I0*Ps + bgd  
    124129         

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/OblateCoreShell_v40.ipf

    r570 r633  
    210210        // local variables 
    211211        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43 
     212        Variable siq,sit 
    212213         
    213214        PI43=4.0/3.0*PI 
    … …  
    220221        vc = PI43*aa*aa*bb 
    221222        vt = PI43*trmaj*trmaj*trmin 
    222         gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc 
    223         gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps 
     223        if(uq == 0) 
     224                siq = 1/3 
     225        else 
     226                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq 
     227        endif 
     228        if(ut == 0) 
     229                sit = 1/3 
     230        else 
     231                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut 
     232        endif 
     233         
     234        gfnc = 3.0*siq*vc*delpc 
     235        gfnt = 3.0*sit*vt*delps 
    224236        tgfn = gfnc+gfnt 
    225237        gfn4 = tgfn*tgfn 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/ProlateCoreShell_v40.ipf

    r570 r633  
    214214        // local variables 
    215215        Variable aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn2,pi43,gfn 
     216        Variable siq,sit 
    216217 
    217218        PI43=4.0/3.0*PI 
    … …  
    224225        vc = PI43*aa*bb*bb 
    225226        vt = PI43*trmaj*trmin*trmin 
    226         gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc 
    227         gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps 
     227         
     228        if(uq == 0.0) 
     229                siq = 1.0/3.0 
     230   else 
     231                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq 
     232        endif 
     233        if(ut == 0.0) 
     234                sit = 1.0/3.0 
     235        else 
     236                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut 
     237        endif 
     238         
     239        gfnc = 3.0*siq*vc*delpc 
     240        gfnt = 3.0*sit*vt*delps 
    228241        gfn = gfnc+gfnt 
    229242        gfn2 = gfn*gfn 

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/StackedDiscs_v40.ipf

    r570 r633  
    221221   //Local variables  
    222222        Variable totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,kk,sqq,dexpt 
     223        Variable be1,be2,si1,si2 
    223224         
    224225        dr1 = rhoc-rhosolv 
    … …  
    233234        sinarg2 = qq*(length+thick)*cos(dum) 
    234235         
    235         t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(bessJ(1,besarg1)/besarg1) 
    236         t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(bessJ(1,besarg2)/besarg2) 
     236        if(besarg1 == 0) 
     237                be1 = 0.5 
     238        else 
     239                be1 = bessJ(1,besarg1)/besarg1 
     240        endif 
     241        if(besarg2 == 0) 
     242                be2 = 0.5 
     243        else 
     244                be2 = bessJ(1,besarg2)/besarg2 
     245        endif 
     246         
     247        if(sinarg1 == 0) 
     248                si1 = 1 
     249        else 
     250                si1 = sin(sinarg1)/sinarg1 
     251        endif 
     252        if(sinarg2 == 0) 
     253                si2 = 1 
     254        else 
     255                si2 = sin(sinarg2)/sinarg2 
     256        endif 
     257         
     258        t1 = 2*area*(2*length)*dr1*si1*be1 
     259        t2 = 2*area*dr2*(totald*si2-2*length*si1)*be2 
    237260         
    238261        retval =((t1+t2)^2)*sin(dum) 
    … …  
    240263        // loop for the structure facture S(q) 
    241264         
    242              kk=1 
    243              sqq=0.0 
    244       do 
     265        kk=1 
     266        sqq=0.0 
     267   do 
    245268                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0 
    246269                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt) 
    247270 
    248                 kk+=1 
     271                kk+=1 
    249272        while (kk<N)                             
    250273         

  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/UniformEllipsoid_v40.ipf

    r570 r633  
    197197        arg = qq*ra*sqrt(1+theta^2*(nu^2-1)) 
    198198         
    199         retval = 9*((sin(arg)-arg*cos(arg))/arg^3)^2 
    200          
     199        if(arg == 0.0) 
     200        retval =1.0/3.0 
     201    else 
     202        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg) 
     203    endif 
     204    retval *= retval 
     205    retval *= 9.0 
     206         
    201207    return retval 
    202208     
Note: See TracChangeset for help on using the changeset viewer.