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.

File:
1 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); 
Note: See TracChangeset for help on using the changeset viewer.