Changeset 975ec8e in sasview for sansmodels/src/libigor/libCylinder.c
- Timestamp:
- Aug 14, 2009 5:58:58 PM (15 years ago)
- 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
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sansmodels/src/libigor/libCylinder.c
r96b59384 r975ec8e 80 80 double va,vb; //upper and lower integration limits 81 81 double summ,zi,yyy,answer,vell; //running tally of integration 82 double summj,vaj,vbj,zij,arg ; //for the inner integration82 double summj,vaj,vbj,zij,arg, si; //for the inner integration 83 83 84 84 Pi = 4.0*atan(1.0); … … 114 114 //now calculate outer integral 115 115 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; 117 122 summ += yyy; 118 123 } … … 155 160 double va,vb; //upper and lower integration limits 156 161 double summ,zi,yyy,answer,vell; //running tally of integration 157 double summj,vaj,vbj,zij,arg ; //for the inner integration162 double summj,vaj,vbj,zij,arg,si; //for the inner integration 158 163 159 164 Pi = 4.0*atan(1.0); … … 189 194 //now calculate outer integral 190 195 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; 192 202 summ += yyy; 193 203 } … … 893 903 return(inten+bkg); 894 904 } 895 905 /* LamellarFFX : calculates the form factor of a lamellar structure - no S(q) effects included 906 -NO polydispersion included 907 */ 908 double 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 } 896 933 /* LamellarPSX : calculates the form factor of a lamellar structure - with S(q) effects included 897 934 ------- … … 1682 1719 { 1683 1720 // 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; 1685 1722 1686 1723 Pi = 4.0*atan(1.0); … … 1695 1732 vc = pi43*aa*bb*bb; 1696 1733 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; 1699 1746 gfn = gfnc+gfnt; 1700 1747 gfn2 = gfn*gfn; … … 1714 1761 { 1715 1762 // 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; 1717 1764 1718 1765 Pi = 4.0*atan(1.0); … … 1726 1773 vc = pi43*aa*aa*bb; 1727 1774 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; 1730 1787 tgfn = gfnc+gfnt; 1731 1788 gfn4 = tgfn*tgfn; … … 1744 1801 1745 1802 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; 1748 1806 vcyl=Pi*radius*radius*zi; 1749 1807 Pq *= vcyl*vcyl; … … 1764 1822 1765 1823 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 } 1767 1827 1768 1828 vcyl=Pi*zi*zi*Lc; … … 1811 1871 // dum is the dummy variable for the integration (theta) 1812 1872 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; 1814 1874 double Pi; 1815 1875 … … 1825 1885 sinarg1 = qq*length*cos(dum); 1826 1886 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; 1830 1910 1831 1911 retval = ((t1+t2)*(t1+t2))*sin(dum); … … 1839 1919 { 1840 1920 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; 1842 1922 double Pi; 1843 1923 … … 1854 1934 sinarg2 = qq*(length+thick)*cos(dum); 1855 1935 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; 1858 1959 1859 1960 retval = ((t1+t2)*(t1+t2))*sin(dum); … … 1906 2007 1907 2008 //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; 1909 2010 double Pi; 1910 2011 int kk; … … 1923 2024 sinarg2 = qq*(length+thick)*cos(dum); 1924 2025 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); 1927 2049 1928 2050 retval =((t1+t2)*(t1+t2))*sin(dum); … … 1937 2059 // end of loop for S(q) 1938 2060 sqq=1.0+2.0*sqq/N; 2061 1939 2062 retval *= sqq; 1940 2063 … … 2016 2139 nu = nua/a; 2017 2140 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 } 2020 2146 retval *= retval; 2021 2147 retval *= 9; … … 2032 2158 arg1 = qq*rshell*sqrt(1-dum*dum); //1=shell (outer radius) 2033 2159 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. 2036 2171 psi = 1/(1-gamma*gamma)*(lam1 - gamma*gamma*lam2); //SRK 10/19/00 2037 2038 2172 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 } 2040 2178 2041 2179 retval = psi*psi*t2*t2; … … 2084 2222 arg += cc*cc*dy*dy; 2085 2223 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 } 2088 2229 return (val); 2089 2230 … … 2099 2240 // h is the HALF-LENGTH of the cylinder = L/2 (A) 2100 2241 2101 double besarg,bj,retval,d1,t1,b1,t2,b2 ; //Local variables2242 double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si; //Local variables 2102 2243 2103 2244 2104 2245 besarg = qq*rr*sin(theta); 2105 2246 siarg = qq * h * cos(theta); 2106 2247 bj =NR_BessJ1(besarg); 2107 2248 2108 2249 //* Computing 2nd power */ 2109 d1 = sin( qq * h * cos(theta));2250 d1 = sin(siarg); 2110 2251 t1 = d1 * d1; 2111 2252 //* Computing 2nd power */ … … 2113 2254 t2 = d1 * d1 * 4.0 * sin(theta); 2114 2255 //* Computing 2nd power */ 2115 d1 = qq * h * cos(theta);2256 d1 = siarg; 2116 2257 b1 = d1 * d1; 2117 2258 //* Computing 2nd power */ 2118 2259 d1 = qq * rr * sin(theta); 2119 2260 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; 2121 2272 2122 2273 return (retval); … … 2136 2287 2137 2288 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 } 2140 2294 2141 2295 //square it
Note: See TracChangeset
for help on using the changeset viewer.