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

Fixed a few more models to return proper numerical values for I(q=0). All that are "fixable" hve been fixed. Ones that are not fixable are either returning NaN or suffer from numerical instabilities at extreme input conditions. Beware.

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

Legend:

Unmodified
Added
Removed

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

    r632 r634  
    19191919                si1 = sin(sinarg1)/sinarg1; 
    19201920        } 
    1921         if (besarg2 == 0.0){ 
     1921        if (sinarg2 == 0.0){ 
    19221922                si2 = 1.0; 
    19231923        }else{ 
    … …  
    19681968                si1 = sin(sinarg1)/sinarg1; 
    19691969        } 
    1970         if (besarg2 == 0.0){ 
     1970        if (sinarg2 == 0.0){ 
    19711971                si2 = 1.0; 
    19721972        }else{ 
    … …  
    20582058                si1 = sin(sinarg1)/sinarg1; 
    20592059        } 
    2060         if (besarg2 == 0.0){ 
     2060        if (sinarg2 == 0.0){ 
    20612061                si2 = 1.0; 
    20622062        }else{ 
    … …  
    24562456        double summ,zi,yyy,answer;                      //running tally of integration 
    24572457        double summj,vaj,vbj,zij;                       //for the inner integration 
    2458         double SphCyl_tmp[7],arg1,arg2,inner; 
     2458        double SphCyl_tmp[7],arg1,arg2,inner,be; 
    24592459         
    24602460         
    … …  
    25062506                arg2 = x*rad*sin(zi); 
    25072507                yyy = inner; 
     2508 
     2509                if(arg2 == 0) { 
     2510                        be = 0.5; 
     2511                } else { 
     2512                        be = NR_BessJ1(arg2)/arg2; 
     2513                } 
    25082514                 
    25092515                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    2510                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
     2516                        yyy += Pi*rad*rad*len*2.0*be; 
    25112517                } else { 
    2512                         yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 
     2518                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 
    25132519                } 
    25142520                yyy *= yyy; 
    … …  
    25372543        double val,arg1,arg2; 
    25382544        double scale,bkg,sldc,slds; 
    2539         double len,rad,hDist,endRad; 
     2545        double len,rad,hDist,endRad,be; 
    25402546        scale = w[0]; 
    25412547        rad = w[1]; 
    … …  
    25512557        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 
    25522558         
    2553         val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; 
     2559        if(arg2 == 0) { 
     2560                be = 0.5; 
     2561        } else { 
     2562                be = NR_BessJ1(arg2)/arg2; 
     2563        } 
     2564        val = cos(arg1)*(1.0-tt*tt)*be; 
    25542565         
    25552566        return(val); 
    … …  
    25732584        double summ,zi,yyy,answer;                      //running tally of integration 
    25742585        double summj,vaj,vbj,zij;                       //for the inner integration 
    2575         double CLens_tmp[7],arg1,arg2,inner,hh; 
     2586        double CLens_tmp[7],arg1,arg2,inner,hh,be; 
    25762587         
    25772588         
    … …  
    26262637                yyy = inner; 
    26272638                 
     2639                if(arg2 == 0) { 
     2640                        be = 0.5; 
     2641                } else { 
     2642                        be = NR_BessJ1(arg2)/arg2; 
     2643                } 
     2644                 
    26282645                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    2629                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
     2646                        yyy += Pi*rad*rad*len*2.0*be; 
    26302647                } else { 
    2631                         yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 
     2648                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 
    26322649                } 
    26332650                yyy *= yyy; 
    … …  
    26672684        double summ,zi,yyy,answer;                      //running tally of integration 
    26682685        double summj,vaj,vbj,zij;                       //for the inner integration 
    2669         double arg1,arg2,inner,hh; 
     2686        double arg1,arg2,inner,hh,be; 
    26702687         
    26712688         
    … …  
    27102727                yyy = inner; 
    27112728                 
     2729                if(arg2 == 0) { 
     2730                        be = 0.5; 
     2731                } else { 
     2732                        be = NR_BessJ1(arg2)/arg2; 
     2733                } 
     2734                 
    27122735                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    2713                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
     2736                        yyy += Pi*rad*rad*len*2.0*be; 
    27142737                } else { 
    2715                         yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 
     2738                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 
    27162739                } 
     2740                 
     2741                 
     2742                 
    27172743                yyy *= yyy; 
    27182744                yyy *= sin(zi);         // = |A(q)|^2*sin(theta) 
    … …  
    27422768        double val,arg1,arg2; 
    27432769        double scale,bkg,sldc,slds; 
    2744         double len,rad,hDist,endRad; 
     2770        double len,rad,hDist,endRad,be; 
    27452771        scale = w[0]; 
    27462772        rad = w[1]; 
    … …  
    27562782        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 
    27572783         
    2758         val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; 
     2784        if(arg2 == 0) { 
     2785                be = 0.5; 
     2786        } else { 
     2787                be = NR_BessJ1(arg2)/arg2; 
     2788        } 
     2789         
     2790        val = cos(arg1)*(1.0-tt*tt)*be; 
    27592791         
    27602792        return(val); 
    … …  
    27782810        double summ,zi,yyy,answer;                      //running tally of integration 
    27792811        double summj,vaj,vbj,zij;                       //for the inner integration 
    2780         double Dumb_tmp[7],arg1,arg2,inner; 
     2812        double Dumb_tmp[7],arg1,arg2,inner,be; 
    27812813         
    27822814         
    … …  
    28312863                yyy = inner; 
    28322864                 
     2865                if(arg2 == 0) { 
     2866                        be = 0.5; 
     2867                } else { 
     2868                        be = NR_BessJ1(arg2)/arg2; 
     2869                } 
     2870                 
    28332871                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    2834                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
     2872                        yyy += Pi*rad*rad*len*2.0*be; 
    28352873                } else { 
    2836                         yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 
     2874                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 
    28372875                } 
    28382876                yyy *= yyy; 
    … …  
    28722910        double summ,zi,yyy,answer;                      //running tally of integration 
    28732911        double summj,vaj,vbj,zij;                       //for the inner integration 
    2874         double arg1,arg2,inner; 
     2912        double arg1,arg2,inner,be; 
    28752913         
    28762914         
    … …  
    29152953                yyy = inner; 
    29162954                 
     2955                if(arg2 == 0) { 
     2956                        be = 0.5; 
     2957                } else { 
     2958                        be = NR_BessJ1(arg2)/arg2; 
     2959                } 
     2960                 
    29172961                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h 
    2918                         yyy += Pi*rad*rad*len*2.0*NR_BessJ1(arg2)/arg2; 
     2962                        yyy += Pi*rad*rad*len*2.0*be; 
    29192963                } else { 
    2920                         yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*NR_BessJ1(arg2)/arg2; 
     2964                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be; 
    29212965                } 
    29222966                yyy *= yyy; 
    … …  
    29482992        double val,arg1,arg2; 
    29492993        double scale,bkg,sldc,slds; 
    2950         double len,rad,hDist,endRad; 
     2994        double len,rad,hDist,endRad,be; 
    29512995        scale = w[0]; 
    29522996        rad = w[1]; 
    … …  
    29623006        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt); 
    29633007         
    2964         val = cos(arg1)*(1.0-tt*tt)*NR_BessJ1(arg2)/arg2; 
     3008        if(arg2 == 0) { 
     3009                be = 0.5; 
     3010        } else { 
     3011                be = NR_BessJ1(arg2)/arg2; 
     3012        } 
     3013        val = cos(arg1)*(1.0-tt*tt)*be; 
    29653014         
    29663015        return(val); 
    … …  
    30593108        double sinarg1,sinarg2; 
    30603109        double t1,t2,t3; 
    3061         double retval; 
     3110        double retval,si1,si2,be1,be2; 
    30623111         
    30633112        double Pi = 4.0*atan(1.0); 
    … …  
    30743123        sinarg2 = qq*(length+facthick)*cos(dum); 
    30753124         
    3076         t1 = 2.0*vol1*dr1*sin(sinarg1)/sinarg1*NR_BessJ1(besarg1)/besarg1; 
    3077         t2 = 2.0*vol2*dr2*sin(sinarg2)/sinarg2*NR_BessJ1(besarg2)/besarg2; 
    3078         t3 = 2.0*vol3*dr3*sin(sinarg2)/sinarg2*NR_BessJ1(besarg1)/besarg1; 
     3125        if(besarg1 == 0) { 
     3126                be1 = 0.5; 
     3127        } else { 
     3128                be1 = NR_BessJ1(besarg1)/besarg1; 
     3129        } 
     3130        if(besarg2 == 0) { 
     3131                be2 = 0.5; 
     3132        } else { 
     3133                be2 = NR_BessJ1(besarg2)/besarg2; 
     3134        }        
     3135        if(sinarg1 == 0) { 
     3136                si1 = 1.0; 
     3137        } else { 
     3138                si1 = sin(sinarg1)/sinarg1; 
     3139        } 
     3140        if(sinarg2 == 0) { 
     3141                si2 = 1.0; 
     3142        } else { 
     3143                si2 = sin(sinarg2)/sinarg2; 
     3144        } 
     3145        t1 = 2.0*vol1*dr1*si1*be1; 
     3146        t2 = 2.0*vol2*dr2*si2*be2; 
     3147        t3 = 2.0*vol3*dr3*si2*be1; 
    30793148         
    30803149        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum); 

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

    r632 r634  
    12871287        qr=x*rcore; 
    12881288        contr = rhocore-rhoshel; 
    1289         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1289        if(qr == 0){ 
     1290                bes = 1.0; 
     1291        }else{ 
     1292                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1293        } 
    12901294        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    12911295        f = vol*bes*contr; 
    … …  
    12931297        qr=x*(rcore+thick); 
    12941298        contr = rhoshel-rhosolv; 
    1295         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1299        if(qr == 0){ 
     1300                bes = 1.0; 
     1301        }else{ 
     1302                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1303        } 
    12961304        vol = 4.0*pi/3.0*pow((rcore+thick),3); 
    12971305        f += vol*bes*contr; 
    … …  
    13421350        qr=x*rcore; 
    13431351        contr = rhocore-rhoshel1; 
    1344         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1352        if(qr == 0){ 
     1353                bes = 1.0; 
     1354        }else{ 
     1355                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1356        } 
    13451357        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    13461358        f = vol*bes*contr; 
    … …  
    13481360        qr=x*(rcore+thick1); 
    13491361        contr = rhoshel1-rhoshel2; 
    1350         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1362        if(qr == 0){ 
     1363                bes = 1.0; 
     1364        }else{ 
     1365                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1366        } 
    13511367        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
    13521368        f += vol*bes*contr; 
    … …  
    13541370        qr=x*(rcore+thick1+thick2); 
    13551371        contr = rhoshel2-rhosolv; 
    1356         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1372        if(qr == 0){ 
     1373                bes = 1.0; 
     1374        }else{ 
     1375                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1376        } 
    13571377        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
    13581378        f += vol*bes*contr; 
    … …  
    14091429        qr=x*rcore; 
    14101430        contr = rhocore-rhoshel1; 
    1411         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1431        if(qr == 0){ 
     1432                bes = 1.0; 
     1433        }else{ 
     1434                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1435        } 
    14121436        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    14131437        f = vol*bes*contr; 
    … …  
    14151439        qr=x*(rcore+thick1); 
    14161440        contr = rhoshel1-rhoshel2; 
    1417         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1441        if(qr == 0){ 
     1442                bes = 1.0; 
     1443        }else{ 
     1444                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1445        } 
    14181446        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
    14191447        f += vol*bes*contr; 
    … …  
    14211449        qr=x*(rcore+thick1+thick2); 
    14221450        contr = rhoshel2-rhoshel3; 
    1423         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1451        if(qr == 0){ 
     1452                bes = 1.0; 
     1453        }else{ 
     1454                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1455        } 
    14241456        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
    14251457        f += vol*bes*contr; 
    … …  
    14271459        qr=x*(rcore+thick1+thick2+thick3); 
    14281460        contr = rhoshel3-rhosolv; 
    1429         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1461        if(qr == 0){ 
     1462                bes = 1.0; 
     1463        }else{ 
     1464                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1465        } 
    14301466        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
    14311467        f += vol*bes*contr; 
    … …  
    14851521        qr=x*rcore; 
    14861522        contr = rhocore-rhoshel1; 
    1487         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1523        if(qr == 0){ 
     1524                bes = 1.0; 
     1525        }else{ 
     1526                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1527        } 
    14881528        vol = 4.0*pi/3.0*rcore*rcore*rcore; 
    14891529        f = vol*bes*contr; 
    … …  
    14911531        qr=x*(rcore+thick1); 
    14921532        contr = rhoshel1-rhoshel2; 
    1493         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1533        if(qr == 0){ 
     1534                bes = 1.0; 
     1535        }else{ 
     1536                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1537        } 
    14941538        vol = 4.0*pi/3.0*(rcore+thick1)*(rcore+thick1)*(rcore+thick1); 
    14951539        f += vol*bes*contr; 
    … …  
    14971541        qr=x*(rcore+thick1+thick2); 
    14981542        contr = rhoshel2-rhoshel3; 
    1499         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1543        if(qr == 0){ 
     1544                bes = 1.0; 
     1545        }else{ 
     1546                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1547        } 
    15001548        vol = 4.0*pi/3.0*(rcore+thick1+thick2)*(rcore+thick1+thick2)*(rcore+thick1+thick2); 
    15011549        f += vol*bes*contr; 
    … …  
    15031551        qr=x*(rcore+thick1+thick2+thick3); 
    15041552        contr = rhoshel3-rhoshel4; 
    1505         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1553        if(qr == 0){ 
     1554                bes = 1.0; 
     1555        }else{ 
     1556                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1557        } 
    15061558        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3)*(rcore+thick1+thick2+thick3); 
    15071559        f += vol*bes*contr; 
    … …  
    15091561        qr=x*(rcore+thick1+thick2+thick3+thick4); 
    15101562        contr = rhoshel4-rhosolv; 
    1511         bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1563        if(qr == 0){ 
     1564                bes = 1.0; 
     1565        }else{ 
     1566                bes = 3.0*(sin(qr)-qr*cos(qr))/(qr*qr*qr); 
     1567        } 
    15121568        vol = 4.0*pi/3.0*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4)*(rcore+thick1+thick2+thick3+thick4); 
    15131569        f += vol*bes*contr; 

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

    r632 r634  
    154154        ans = G1*exp(-x*x*Rg1*Rg1/3.0); 
    155155        ans += B1*pow((erf1*erf1*erf1/x),Pow1); 
     156         
     157        if(x == 0) { 
     158                ans = G1; 
     159        } 
    156160         
    157161        ans *= scale; 
    … …  
    190194        ans += G2*exp(-x*x*Rg2*Rg2/3.0); 
    191195        ans += B2*pow((erf2*erf2*erf2/x),Pow2); 
     196 
     197        if(x == 0) { 
     198                ans = G1+G2; 
     199        } 
    192200         
    193201        ans *= scale; 
    … …  
    232240        ans += G2*exp(-x*x*Rg2*Rg2/3.0) + B2*exp(-x*x*Rg3*Rg3/3.0)*pow((erf2*erf2*erf2/x),Pow2); 
    233241        ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*pow((erf3*erf3*erf3/x),Pow3); 
     242 
     243        if(x == 0) { 
     244                ans = G1+G2+G3; 
     245        } 
    234246         
    235247        ans *= scale; 
    … …  
    280292        ans += G3*exp(-x*x*Rg3*Rg3/3.0) + B3*exp(-x*x*Rg4*Rg4/3.0)*pow((erf3*erf3*erf3/x),Pow3); 
    281293        ans += G4*exp(-x*x*Rg4*Rg4/3.0) + B4*pow((erf4*erf4*erf4/x),Pow4); 
     294 
     295        if(x == 0) { 
     296                ans = G1+G2+G3+G4; 
     297        } 
    282298         
    283299        ans *= scale; 
Note: See TracChangeset for help on using the changeset viewer.