Ignore:
Timestamp:
Aug 14, 2009 5:58:58 PM (15 years ago)
Author:
Jae Cho <jhjcho@…>
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:
d5bd424
Parents:
cad821b
Message:

working on 2D models. Still need smore corrections and unit tests.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/libigor/libCylinder.c

    r96b59384 r975ec8e  
    8080        double va,vb;           //upper and lower integration limits 
    8181        double summ,zi,yyy,answer,vell;                 //running tally of integration 
    82         double summj,vaj,vbj,zij,arg;                   //for the inner integration 
     82        double summj,vaj,vbj,zij,arg, si;                       //for the inner integration 
    8383 
    8484        Pi = 4.0*atan(1.0); 
    … …  
    114114                //now calculate outer integral 
    115115                arg = q*length*zi/2; 
    116                 yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg; 
     116                if (arg == 0){ 
     117                        si = 1; 
     118                }else{ 
     119                        si = sin(arg) * sin(arg) / arg / arg; 
     120                } 
     121                yyy = Gauss76Wt[i] * answer * si; 
    117122                summ += yyy; 
    118123        } 
    … …  
    155160        double va,vb;           //upper and lower integration limits 
    156161        double summ,zi,yyy,answer,vell;                 //running tally of integration 
    157         double summj,vaj,vbj,zij,arg;                   //for the inner integration 
     162        double summj,vaj,vbj,zij,arg,si;                        //for the inner integration 
    158163 
    159164        Pi = 4.0*atan(1.0); 
    … …  
    189194                //now calculate outer integral 
    190195                arg = q*length*zi/2; 
    191                 yyy = Gauss76Wt[i] * answer * sin(arg) * sin(arg) / arg / arg; 
     196                if (arg == 0){ 
     197                        si = 1; 
     198                }else{ 
     199                        si = sin(arg) * sin(arg) / arg / arg; 
     200                } 
     201                yyy = Gauss76Wt[i] * answer * si; 
    192202                summ += yyy; 
    193203        } 
    … …  
    893903        return(inten+bkg); 
    894904} 
    895  
     905/*      LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included 
     906                                                -NO polydispersion included 
     907*/ 
     908double lamellar_kernel(double dp[], double q){ 
     909        double scale,del,sld_bi,sld_sol,contr,bkg;              //local variables of coefficient wave 
     910        double inten, qval,Pq; 
     911        double Pi; 
     912 
     913 
     914        Pi = 4.0*atan(1.0); 
     915        scale = dp[0]; 
     916        del = dp[1]; 
     917        sld_bi = dp[2]; 
     918        sld_sol = dp[3]; 
     919        bkg = dp[4]; 
     920        qval = q; 
     921        contr = sld_bi -sld_sol; 
     922 
     923        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)); 
     924 
     925        inten = 2.0*Pi*scale*Pq/(qval*qval);            //this is now dimensionless... 
     926 
     927        inten /= del;                   //normalize by the thickness (in A) 
     928 
     929        inten *= 1.0e8;         // 1/A to 1/cm 
     930 
     931        return(inten+bkg); 
     932} 
    896933/*      LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included 
    897934------- 
    … …  
    16821719{ 
    16831720        // local variables 
    1684         double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,gfn2,pi43,gfn,Pi; 
     1721        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi; 
    16851722 
    16861723        Pi = 4.0*atan(1.0); 
    … …  
    16951732        vc = pi43*aa*bb*bb; 
    16961733        vt = pi43*trmaj*trmin*trmin; 
    1697         gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc; 
    1698         gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps; 
     1734        if (uq == 0){ 
     1735                siq = 1/3; 
     1736        }else{ 
     1737                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; 
     1738        } 
     1739        if (ut == 0){ 
     1740                sit = 1/3; 
     1741        }else{ 
     1742                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; 
     1743        } 
     1744        gfnc = 3.0*siq*vc*delpc; 
     1745        gfnt = 3.0*sit*vt*delps; 
    16991746        gfn = gfnc+gfnt; 
    17001747        gfn2 = gfn*gfn; 
    … …  
    17141761{ 
    17151762        // local variables 
    1716         double aa,bb,u2,ut2,uq,ut,vc,vt,gfnc,gfnt,tgfn,gfn4,pi43,Pi; 
     1763        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi; 
    17171764 
    17181765        Pi = 4.0*atan(1.0); 
    … …  
    17261773        vc = pi43*aa*aa*bb; 
    17271774        vt = pi43*trmaj*trmaj*trmin; 
    1728         gfnc = 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq*vc*delpc; 
    1729         gfnt = 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut*vt*delps; 
     1775        if (uq == 0){ 
     1776                siq = 1/3; 
     1777        }else{ 
     1778                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq; 
     1779        } 
     1780        if (ut == 0){ 
     1781                sit = 1/3; 
     1782        }else{ 
     1783                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut; 
     1784        } 
     1785        gfnc = 3.0*siq*vc*delpc; 
     1786        gfnt = 3.0*sit*vt*delps; 
    17301787        tgfn = gfnc+gfnt; 
    17311788        gfn4 = tgfn*tgfn; 
    … …  
    17441801 
    17451802    Pq = Sk_WR(q,zi,lb);                //does not have cross section term 
    1746     Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
    1747  
     1803    if (qr !=0){ 
     1804        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
     1805    } //else Pk *=1; 
    17481806    vcyl=Pi*radius*radius*zi; 
    17491807    Pq *= vcyl*vcyl; 
    … …  
    17641822 
    17651823    Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term 
    1766     Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
     1824    if (qr !=0){ 
     1825        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr); 
     1826    } 
    17671827 
    17681828    vcyl=Pi*zi*zi*Lc; 
    … …  
    18111871        // dum is the dummy variable for the integration (theta) 
    18121872 
    1813         double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval; 
     1873        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; 
    18141874        double Pi; 
    18151875 
    … …  
    18251885        sinarg1 = qq*length*cos(dum); 
    18261886        sinarg2 = qq*(length+facthick)*cos(dum); 
    1827  
    1828         t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
    1829         t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
     1887        if (besarg1 == 0){ 
     1888                be1 = 0.5; 
     1889        }else{ 
     1890                be1 = NR_BessJ1(besarg1)/besarg1; 
     1891        } 
     1892        if (besarg2 == 0){ 
     1893                be2 = 0.5; 
     1894        }else{ 
     1895                be2 = NR_BessJ1(besarg2)/besarg2; 
     1896        } 
     1897        if (sinarg1 == 0){ 
     1898                si1 = 1; 
     1899        }else{ 
     1900                si1 = sin(sinarg1)/sinarg1; 
     1901        } 
     1902        if (besarg2 == 0){ 
     1903                si2 = 1; 
     1904        }else{ 
     1905                si2 = sin(sinarg2)/sinarg2; 
     1906        } 
     1907 
     1908        t1 = 2.0*vol1*dr1*si1*be1; 
     1909        t2 = 2.0*vol2*dr2*si2*be2; 
    18301910 
    18311911        retval = ((t1+t2)*(t1+t2))*sin(dum); 
    … …  
    18391919{ 
    18401920 
    1841     double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,t1,t2,retval; 
     1921    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval; 
    18421922    double Pi; 
    18431923 
    … …  
    18541934    sinarg2 = qq*(length+thick)*cos(dum); 
    18551935 
    1856     t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
    1857     t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
     1936        if (besarg1 == 0){ 
     1937                be1 = 0.5; 
     1938        }else{ 
     1939                be1 = NR_BessJ1(besarg1)/besarg1; 
     1940        } 
     1941        if (besarg2 == 0){ 
     1942                be2 = 0.5; 
     1943        }else{ 
     1944                be2 = NR_BessJ1(besarg2)/besarg2; 
     1945        } 
     1946        if (sinarg1 == 0){ 
     1947                si1 = 1; 
     1948        }else{ 
     1949                si1 = sin(sinarg1)/sinarg1; 
     1950        } 
     1951        if (besarg2 == 0){ 
     1952                si2 = 1; 
     1953        }else{ 
     1954                si2 = sin(sinarg2)/sinarg2; 
     1955        } 
     1956 
     1957    t1 = 2.0*vol1*dr1*si1*be1; 
     1958    t2 = 2.0*vol2*dr2*si2*be2; 
    18581959 
    18591960    retval = ((t1+t2)*(t1+t2))*sin(dum); 
    … …  
    19062007 
    19072008        //Local variables 
    1908         double totald,dr1,dr2,besarg1,besarg2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; 
     2009        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt; 
    19092010        double Pi; 
    19102011        int kk; 
    … …  
    19232024        sinarg2 = qq*(length+thick)*cos(dum); 
    19242025 
    1925         t1 = 2*area*(2*length)*dr1*(sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg1)/besarg1); 
    1926         t2 = 2*area*dr2*(totald*sin(sinarg2)/sinarg2-2*length*sin(sinarg1)/sinarg1)*(NR_BessJ1(besarg2)/besarg2); 
     2026        if (besarg1 == 0){ 
     2027                be1 = 0.5; 
     2028        }else{ 
     2029                be1 = NR_BessJ1(besarg1)/besarg1; 
     2030        } 
     2031        if (besarg2 == 0){ 
     2032                be2 = 0.5; 
     2033        }else{ 
     2034                be2 = NR_BessJ1(besarg2)/besarg2; 
     2035        } 
     2036        if (sinarg1 == 0){ 
     2037                si1 = 1; 
     2038        }else{ 
     2039                si1 = sin(sinarg1)/sinarg1; 
     2040        } 
     2041        if (besarg2 == 0){ 
     2042                si2 = 1; 
     2043        }else{ 
     2044                si2 = sin(sinarg2)/sinarg2; 
     2045        } 
     2046 
     2047        t1 = 2*area*(2*length)*dr1*(si1)*(be1); 
     2048        t2 = 2*area*dr2*(totald*si2-2*length*si1)*(be2); 
    19272049 
    19282050        retval =((t1+t2)*(t1+t2))*sin(dum); 
    … …  
    19372059        // end of loop for S(q) 
    19382060        sqq=1.0+2.0*sqq/N; 
     2061 
    19392062        retval *= sqq; 
    19402063 
    … …  
    20162139    nu = nua/a; 
    20172140    arg = qq*a*sqrt(1+dum*dum*(nu*nu-1)); 
    2018  
    2019     retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); 
     2141    if (arg == 0){ 
     2142        retval =1/3; 
     2143    }else{ 
     2144        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg); 
     2145    } 
    20202146    retval *= retval; 
    20212147    retval *= 9; 
    … …  
    20322158    arg1 = qq*rshell*sqrt(1-dum*dum);           //1=shell (outer radius) 
    20332159    arg2 = qq*rcore*sqrt(1-dum*dum);                    //2=core (inner radius) 
    2034     lam1 = 2*NR_BessJ1(arg1)/arg1; 
    2035     lam2 = 2*NR_BessJ1(arg2)/arg2; 
     2160    if (arg1 == 0){ 
     2161        lam1 = 1; 
     2162    }else{ 
     2163        lam1 = 2*NR_BessJ1(arg1)/arg1; 
     2164    } 
     2165    if (arg2 == 0){ 
     2166        lam2 = 1; 
     2167    }else{ 
     2168        lam2 = 2*NR_BessJ1(arg2)/arg2; 
     2169    } 
     2170    //Todo: Need to check psi behavior as gamma goes to 1. 
    20362171    psi = 1/(1-gamma*gamma)*(lam1 -  gamma*gamma*lam2);         //SRK 10/19/00 
    2037  
    20382172    sinarg = qq*length*dum/2; 
    2039     t2 = sin(sinarg)/sinarg; 
     2173    if (sinarg == 0){ 
     2174        t2 = 1; 
     2175    }else{ 
     2176        t2 = sin(sinarg)/sinarg; 
     2177    } 
    20402178 
    20412179    retval = psi*psi*t2*t2; 
    … …  
    20842222    arg += cc*cc*dy*dy; 
    20852223    arg = q*sqrt(arg); 
    2086     val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); 
    2087  
     2224    if (arg == 0){ 
     2225        val = 1; // as arg --> 0, val should go to 1.0 
     2226    }else{ 
     2227        val = 9 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ); 
     2228    } 
    20882229    return (val); 
    20892230 
    … …  
    20992240        // h is the HALF-LENGTH of the cylinder = L/2 (A) 
    21002241 
    2101     double besarg,bj,retval,d1,t1,b1,t2,b2;              //Local variables 
     2242    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;          //Local variables 
    21022243 
    21032244 
    21042245    besarg = qq*rr*sin(theta); 
    2105  
     2246    siarg = qq * h * cos(theta); 
    21062247    bj =NR_BessJ1(besarg); 
    21072248 
    21082249        //* Computing 2nd power */ 
    2109     d1 = sin(qq * h * cos(theta)); 
     2250    d1 = sin(siarg); 
    21102251    t1 = d1 * d1; 
    21112252        //* Computing 2nd power */ 
    … …  
    21132254    t2 = d1 * d1 * 4.0 * sin(theta); 
    21142255        //* Computing 2nd power */ 
    2115     d1 = qq * h * cos(theta); 
     2256    d1 = siarg; 
    21162257    b1 = d1 * d1; 
    21172258        //* Computing 2nd power */ 
    21182259    d1 = qq * rr * sin(theta); 
    21192260    b2 = d1 * d1; 
    2120     retval = t1 * t2 / b1 / b2; 
     2261    if (besarg == 0){ 
     2262        be = sin(theta); 
     2263    }else{ 
     2264        be = t2 / b2; 
     2265    } 
     2266    if (siarg == 0){ 
     2267        si = 1; 
     2268    }else{ 
     2269        si = t1 / b1; 
     2270    } 
     2271    retval = be * si; 
    21212272 
    21222273    return (retval); 
    … …  
    21362287 
    21372288        arg = qq*ra*sqrt((1+nu*nu)/2+(1-nu*nu)*cos(theta)/2); 
    2138  
    2139         retval = 2*NR_BessJ1(arg)/arg; 
     2289        if (arg == 0){ 
     2290                retval = 1; 
     2291        }else{ 
     2292                retval = 2*NR_BessJ1(arg)/arg; 
     2293        } 
    21402294 
    21412295        //square it 
Note: See TracChangeset for help on using the changeset viewer.