Changeset 5f0dcab in sasview for sansmodels/src


Ignore:
Timestamp:
Dec 7, 2011 3:39:05 PM (13 years ago)
Author:
Mathieu Doucet <doucetm@…>
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:
feadd6f
Parents:
2605da22
Message:

Re #4 Get rid of warnings by putting static and defs in the right place.

Location:
sansmodels/src/sans/models/libigor
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/sans/models/libigor/libCylinder.c

    r34c2649 r5f0dcab  
    88#include "GaussWeights.h" 
    99#include "libCylinder.h" 
     10 
     11/////////functions for WRC implementation of flexible cylinders 
     12static double 
     13gammaln(double xx) 
     14{ 
     15  double x,y,tmp,ser; 
     16  static double cof[6]={76.18009172947146,-86.50532032941677, 
     17    24.01409824083091,-1.231739572450155, 
     18    0.1208650973866179e-2,-0.5395239384953e-5}; 
     19  int j; 
     20 
     21  y=x=xx; 
     22  tmp=x+5.5; 
     23  tmp -= (x+0.5)*log(tmp); 
     24  ser=1.000000000190015; 
     25  for (j=0;j<=5;j++) ser += cof[j]/++y; 
     26  return -tmp+log(2.5066282746310005*ser/x); 
     27} 
     28/////////functions for WRC implementation of flexible cylinders 
     29 
     30// 
     31static double 
     32AlphaSquare(double x) 
     33{ 
     34    double yy; 
     35    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
     36 
     37    return (yy); 
     38} 
     39 
     40// 
     41static double 
     42Rgsquarezero(double q, double L, double b) 
     43{ 
     44    double yy; 
     45    yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b)))); 
     46 
     47    return (yy); 
     48} 
     49 
     50// 
     51static double 
     52Rgsquareshort(double q, double L, double b) 
     53{ 
     54    double yy; 
     55    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
     56 
     57    return (yy); 
     58} 
     59 
     60// 
     61static double 
     62Rgsquare(double q, double L, double b) 
     63{ 
     64    double yy; 
     65    yy = AlphaSquare(L/b)*L*b/6.0; 
     66 
     67    return (yy); 
     68} 
     69 
     70// ?? funciton is not used - but should the log actually be log10??? 
     71static double 
     72miu(double x) 
     73{ 
     74    double yy; 
     75    yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x))); 
     76 
     77    return (yy); 
     78} 
     79 
     80//WR named this w (too generic) 
     81static double 
     82w_WR(double x) 
     83{ 
     84    double yy; 
     85    yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 
     86 
     87    return (yy); 
     88} 
     89 
     90// 
     91static double 
     92u1(double q, double L, double b) 
     93{ 
     94    double yy; 
     95 
     96    yy = Rgsquareshort(q,L,b)*q*q; 
     97 
     98    return (yy); 
     99} 
     100 
     101// was named u 
     102static double 
     103u_WR(double q, double L, double b) 
     104{ 
     105    double yy; 
     106    yy = Rgsquare(q,L,b)*q*q; 
     107    return (yy); 
     108} 
     109 
     110 
     111// 
     112static double 
     113Sdebye1(double q, double L, double b) 
     114{ 
     115    double yy; 
     116    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 
     117 
     118    return (yy); 
     119} 
     120 
     121 
     122// 
     123static double 
     124Sdebye(double q, double L, double b) 
     125{ 
     126    double yy; 
     127    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 
     128 
     129    return (yy); 
     130} 
     131 
     132// 
     133static double 
     134Sexv(double q, double L, double b) 
     135{ 
     136    double yy,C1,C2,C3,miu,Rg2; 
     137    C1=1.22; 
     138    C2=0.4288; 
     139    C3=-1.651; 
     140    miu = 0.585; 
     141 
     142    Rg2 = Rgsquare(q,L,b); 
     143 
     144    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 
     145 
     146    return (yy); 
     147} 
     148 
     149// this must be WR modified version 
     150static double 
     151Sexvnew(double q, double L, double b) 
     152{ 
     153    double yy,C1,C2,C3,miu; 
     154    double del=1.05,C_star2,Rg2; 
     155 
     156    C1=1.22; 
     157    C2=0.4288; 
     158    C3=-1.651; 
     159    miu = 0.585; 
     160 
     161    //calculating the derivative to decide on the corection (cutoff) term? 
     162    // I have modified this from WRs original code 
     163 
     164    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 
     165        C_star2 = 0.0; 
     166    } else { 
     167        C_star2 = 1.0; 
     168    } 
     169 
     170    Rg2 = Rgsquare(q,L,b); 
     171 
     172  yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 
     173 
     174    return (yy); 
     175} 
     176 
     177// these are the messy ones 
     178static double 
     179a2short(double q, double L, double b, double p1short, double p2short, double q0) 
     180{ 
     181    double yy,Rg2_sh; 
     182    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 
     183    double pi; 
     184 
     185    E = 2.718281828459045091; 
     186  pi = 4.0*atan(1.0); 
     187    Rg2_sh = Rgsquareshort(q,L,b); 
     188    Rg2_sh2 = Rg2_sh*Rg2_sh; 
     189    t1 = ((q0*q0*Rg2_sh)/(b*b)); 
     190    Et1 = pow(E,t1); 
     191    Emt1 =pow(E,-t1); 
     192    q02 = q0*q0; 
     193    q0p = pow(q0,(-4.0 + p2short) ); 
     194 
     195    //E is the number e 
     196    yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2))))))); 
     197 
     198    return (yy); 
     199} 
     200 
     201// 
     202static double 
     203a1short(double q, double L, double b, double p1short, double p2short, double q0) 
     204{ 
     205    double yy,Rg2_sh; 
     206    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 
     207  double pi; 
     208 
     209    E = 2.718281828459045091; 
     210  pi = 4.0*atan(1.0); 
     211    Rg2_sh = Rgsquareshort(q,L,b); 
     212    Rg2_sh2 = Rg2_sh*Rg2_sh; 
     213    b3 = b*b*b; 
     214    t1 = ((q0*q0*Rg2_sh)/(b*b)); 
     215    Et1 = pow(E,t1); 
     216    Emt1 =pow(E,-t1); 
     217    q02 = q0*q0; 
     218    q0p = pow(q0,(-4.0 + p1short) ); 
     219 
     220    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); 
     221 
     222    return(yy); 
     223} 
     224 
     225 
     226//need to define this on my own 
     227static double 
     228sech_WR(double x) 
     229{ 
     230  return(1/cosh(x)); 
     231} 
     232 
     233// this one will be lots of trouble 
     234static double 
     235a2long(double q, double L, double b, double p1, double p2, double q0) 
     236{ 
     237    double yy,C1,C2,C3,C4,C5,miu,C,Rg2; 
     238    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 
     239    double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 
     240 
     241    pi = 4.0*atan(1.0); 
     242    E = 2.718281828459045091; 
     243    if( L/b > 10.0) { 
     244        C = 3.06/pow((L/b),0.44); 
     245    } else { 
     246        C = 1.0; 
     247    } 
     248 
     249    C1 = 1.22; 
     250    C2 = 0.4288; 
     251    C3 = -1.651; 
     252    C4 = 1.523; 
     253    C5 = 0.1477; 
     254    miu = 0.585; 
     255 
     256    Rg2 = Rgsquare(q,L,b); 
     257    Rg22 = Rg2*Rg2; 
     258    b2 = b*b; 
     259    b3 = b*b*b; 
     260    b4 = b3*b; 
     261    q02 = q0*q0; 
     262    q03 = q0*q0*q0; 
     263    q04 = q03*q0; 
     264    q05 = q04*q0; 
     265 
     266    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
     267 
     268    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; 
     269 
     270    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5); 
     271 
     272    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 
     273 
     274    t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
     275 
     276    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 
     277 
     278    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 
     279 
     280    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     281 
     282    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; 
     283 
     284    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
     285 
     286 
     287    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 
     288 
     289    return (yy); 
     290} 
     291// 
     292static double 
     293a1long(double q, double L, double b, double p1, double p2, double q0) 
     294{ 
     295    double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 
     296    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 
     297    double E,pi; 
     298    double b2,b3,b4,q02,q03,q04,q05,Rg22; 
     299 
     300    pi = 4.0*atan(1.0); 
     301    E = 2.718281828459045091; 
     302 
     303    if( L/b > 10.0) { 
     304        C = 3.06/pow((L/b),0.44); 
     305    } else { 
     306        C = 1.0; 
     307    } 
     308 
     309    C1 = 1.22; 
     310    C2 = 0.4288; 
     311    C3 = -1.651; 
     312    C4 = 1.523; 
     313    C5 = 0.1477; 
     314    miu = 0.585; 
     315 
     316    Rg2 = Rgsquare(q,L,b); 
     317    Rg22 = Rg2*Rg2; 
     318    b2 = b*b; 
     319    b3 = b*b*b; 
     320    b4 = b3*b; 
     321    q02 = q0*q0; 
     322    q03 = q0*q0*q0; 
     323    q04 = q03*q0; 
     324    q05 = q04*q0; 
     325 
     326    t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 
     327 
     328    t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
     329 
     330    t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 
     331 
     332    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     333 
     334    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 
     335 
     336    t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b))); 
     337 
     338    t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
     339 
     340    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
     341 
     342    t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
     343 
     344    t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
     345 
     346    t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 
     347 
     348    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
     349 
     350    t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 
     351 
     352    t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
     353 
     354    t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 
     355 
     356 
     357    yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))); 
     358 
     359    return (yy); 
     360} 
     361 
     362 
     363static double 
     364Sk_WR(double q, double L, double b) 
     365{ 
     366  // 
     367  double p1,p2,p1short,p2short,q0; 
     368  double C,ans,q0short,Sexvmodify,pi; 
     369 
     370  pi = 4.0*atan(1.0); 
     371 
     372  p1 = 4.12; 
     373  p2 = 4.42; 
     374  p1short = 5.36; 
     375  p2short = 5.62; 
     376  q0 = 3.1; 
     377  // 
     378  q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
     379 
     380  // 
     381  if(L/b > 10.0) { 
     382    C = 3.06/pow((L/b),0.44); 
     383  } else { 
     384    C = 1.0; 
     385  } 
     386  // 
     387 
     388  if( L > 4*b ) { // Longer Chains 
     389    if (q*b <= 3.1) {   //Modified by Yun on Oct. 15, 
     390      Sexvmodify = Sexvnew(q, L, b); 
     391      ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
     392    } else { //q(i)*b > 3.1 
     393      ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L); 
     394    } 
     395  } else { //L <= 4*b Shorter Chains 
     396    if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
     397      if (q*b<=0.01) { 
     398        ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
     399      } else { 
     400        ans = Sdebye1(q,L,b); 
     401      } 
     402    } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
     403      ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); 
     404    } 
     405  } 
     406 
     407  return(ans); 
     408  //return(a2long(q, L, b, p1, p2, q0)); 
     409} 
    10410 
    11411/*      CylinderForm  :  calculates the form factor of a cylinder at the give x-value p->x 
    … …  
    20420CylinderForm(double dp[], double q) 
    21421{ 
     422 
    22423        int i; 
    23424        double Pi; 

  • sansmodels/src/sans/models/libigor/libCylinder.h

    r6a2c931 r5f0dcab  
    6868double CSPPKernel(double dp[], double mu, double uu); 
    6969 
    70 /////////functions for WRC implementation of flexible cylinders 
    71 static double 
    72 gammaln(double xx) 
    73 { 
    74         double x,y,tmp,ser; 
    75         static double cof[6]={76.18009172947146,-86.50532032941677, 
    76                 24.01409824083091,-1.231739572450155, 
    77                 0.1208650973866179e-2,-0.5395239384953e-5}; 
    78         int j; 
    7970 
    80         y=x=xx; 
    81         tmp=x+5.5; 
    82         tmp -= (x+0.5)*log(tmp); 
    83         ser=1.000000000190015; 
    84         for (j=0;j<=5;j++) ser += cof[j]/++y; 
    85         return -tmp+log(2.5066282746310005*ser/x); 
    86 } 
    87  
    88 // 
    89 static double 
    90 AlphaSquare(double x) 
    91 { 
    92     double yy; 
    93     yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) ); 
    94  
    95     return (yy); 
    96 } 
    97  
    98 // 
    99 static double 
    100 Rgsquarezero(double q, double L, double b) 
    101 { 
    102     double yy; 
    103     yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b)))); 
    104  
    105     return (yy); 
    106 } 
    107  
    108 // 
    109 static double 
    110 Rgsquareshort(double q, double L, double b) 
    111 { 
    112     double yy; 
    113     yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b); 
    114  
    115     return (yy); 
    116 } 
    117  
    118 // 
    119 static double 
    120 Rgsquare(double q, double L, double b) 
    121 { 
    122     double yy; 
    123     yy = AlphaSquare(L/b)*L*b/6.0; 
    124  
    125     return (yy); 
    126 } 
    127  
    128 // ?? funciton is not used - but should the log actually be log10??? 
    129 static double 
    130 miu(double x) 
    131 { 
    132     double yy; 
    133     yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x))); 
    134  
    135     return (yy); 
    136 } 
    137  
    138 //WR named this w (too generic) 
    139 static double 
    140 w_WR(double x) 
    141 { 
    142     double yy; 
    143     yy = 0.5*(1 + tanh((x - 1.523)/0.1477)); 
    144  
    145     return (yy); 
    146 } 
    147  
    148 // 
    149 static double 
    150 u1(double q, double L, double b) 
    151 { 
    152     double yy; 
    153  
    154     yy = Rgsquareshort(q,L,b)*q*q; 
    155  
    156     return (yy); 
    157 } 
    158  
    159 // was named u 
    160 static double 
    161 u_WR(double q, double L, double b) 
    162 { 
    163     double yy; 
    164     yy = Rgsquare(q,L,b)*q*q; 
    165     return (yy); 
    166 } 
    167  
    168  
    169 // 
    170 static double 
    171 Sdebye1(double q, double L, double b) 
    172 { 
    173     double yy; 
    174     yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) ); 
    175  
    176     return (yy); 
    177 } 
    178  
    179  
    180 // 
    181 static double 
    182 Sdebye(double q, double L, double b) 
    183 { 
    184     double yy; 
    185     yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2)); 
    186  
    187     return (yy); 
    188 } 
    189  
    190 // 
    191 static double 
    192 Sexv(double q, double L, double b) 
    193 { 
    194     double yy,C1,C2,C3,miu,Rg2; 
    195     C1=1.22; 
    196     C2=0.4288; 
    197     C3=-1.651; 
    198     miu = 0.585; 
    199  
    200     Rg2 = Rgsquare(q,L,b); 
    201  
    202     yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 
    203  
    204     return (yy); 
    205 } 
    206  
    207 // this must be WR modified version 
    208 static double 
    209 Sexvnew(double q, double L, double b) 
    210 { 
    211     double yy,C1,C2,C3,miu; 
    212     double del=1.05,C_star2,Rg2; 
    213  
    214     C1=1.22; 
    215     C2=0.4288; 
    216     C3=-1.651; 
    217     miu = 0.585; 
    218  
    219     //calculating the derivative to decide on the corection (cutoff) term? 
    220     // I have modified this from WRs original code 
    221  
    222     if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) { 
    223         C_star2 = 0.0; 
    224     } else { 
    225         C_star2 = 1.0; 
    226     } 
    227  
    228     Rg2 = Rgsquare(q,L,b); 
    229  
    230         yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu))); 
    231  
    232     return (yy); 
    233 } 
    234  
    235 // these are the messy ones 
    236 static double 
    237 a2short(double q, double L, double b, double p1short, double p2short, double q0) 
    238 { 
    239     double yy,Rg2_sh; 
    240     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p; 
    241     double pi; 
    242  
    243     E = 2.718281828459045091; 
    244         pi = 4.0*atan(1.0); 
    245     Rg2_sh = Rgsquareshort(q,L,b); 
    246     Rg2_sh2 = Rg2_sh*Rg2_sh; 
    247     t1 = ((q0*q0*Rg2_sh)/(b*b)); 
    248     Et1 = pow(E,t1); 
    249     Emt1 =pow(E,-t1); 
    250     q02 = q0*q0; 
    251     q0p = pow(q0,(-4.0 + p2short) ); 
    252  
    253     //E is the number e 
    254     yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2))))))); 
    255  
    256     return (yy); 
    257 } 
    258  
    259 // 
    260 static double 
    261 a1short(double q, double L, double b, double p1short, double p2short, double q0) 
    262 { 
    263     double yy,Rg2_sh; 
    264     double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3; 
    265         double pi; 
    266  
    267     E = 2.718281828459045091; 
    268         pi = 4.0*atan(1.0); 
    269     Rg2_sh = Rgsquareshort(q,L,b); 
    270     Rg2_sh2 = Rg2_sh*Rg2_sh; 
    271     b3 = b*b*b; 
    272     t1 = ((q0*q0*Rg2_sh)/(b*b)); 
    273     Et1 = pow(E,t1); 
    274     Emt1 =pow(E,-t1); 
    275     q02 = q0*q0; 
    276     q0p = pow(q0,(-4.0 + p1short) ); 
    277  
    278     yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2)))))); 
    279  
    280     return(yy); 
    281 } 
    282  
    283  
    284 //need to define this on my own 
    285 static double 
    286 sech_WR(double x) 
    287 { 
    288         return(1/cosh(x)); 
    289 } 
    290  
    291 // this one will be lots of trouble 
    292 static double 
    293 a2long(double q, double L, double b, double p1, double p2, double q0) 
    294 { 
    295     double yy,C1,C2,C3,C4,C5,miu,C,Rg2; 
    296     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi; 
    297     double E,b2,b3,b4,q02,q03,q04,q05,Rg22; 
    298  
    299     pi = 4.0*atan(1.0); 
    300     E = 2.718281828459045091; 
    301     if( L/b > 10.0) { 
    302         C = 3.06/pow((L/b),0.44); 
    303     } else { 
    304         C = 1.0; 
    305     } 
    306  
    307     C1 = 1.22; 
    308     C2 = 0.4288; 
    309     C3 = -1.651; 
    310     C4 = 1.523; 
    311     C5 = 0.1477; 
    312     miu = 0.585; 
    313  
    314     Rg2 = Rgsquare(q,L,b); 
    315     Rg22 = Rg2*Rg2; 
    316     b2 = b*b; 
    317     b3 = b*b*b; 
    318     b4 = b3*b; 
    319     q02 = q0*q0; 
    320     q03 = q0*q0*q0; 
    321     q04 = q03*q0; 
    322     q05 = q04*q0; 
    323  
    324     t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) )); 
    325  
    326     t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L; 
    327  
    328     t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5); 
    329  
    330     t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22); 
    331  
    332     t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
    333  
    334     t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22); 
    335  
    336     t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 
    337  
    338     t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    339  
    340     t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L; 
    341  
    342     t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22); 
    343  
    344  
    345     yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))); 
    346  
    347     return (yy); 
    348 } 
    349 // 
    350 static double 
    351 a1long(double q, double L, double b, double p1, double p2, double q0) 
    352 { 
    353     double yy,C,C1,C2,C3,C4,C5,miu,Rg2; 
    354     double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15; 
    355     double E,pi; 
    356     double b2,b3,b4,q02,q03,q04,q05,Rg22; 
    357  
    358     pi = 4.0*atan(1.0); 
    359     E = 2.718281828459045091; 
    360  
    361     if( L/b > 10.0) { 
    362         C = 3.06/pow((L/b),0.44); 
    363     } else { 
    364         C = 1.0; 
    365     } 
    366  
    367     C1 = 1.22; 
    368     C2 = 0.4288; 
    369     C3 = -1.651; 
    370     C4 = 1.523; 
    371     C5 = 0.1477; 
    372     miu = 0.585; 
    373  
    374     Rg2 = Rgsquare(q,L,b); 
    375     Rg22 = Rg2*Rg2; 
    376     b2 = b*b; 
    377     b3 = b*b*b; 
    378     b4 = b3*b; 
    379     q02 = q0*q0; 
    380     q03 = q0*q0*q0; 
    381     q04 = q03*q0; 
    382     q05 = q04*q0; 
    383  
    384     t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 
    385  
    386     t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    387  
    388     t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 
    389  
    390     t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    391  
    392     t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)))); 
    393  
    394     t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b))); 
    395  
    396     t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
    397  
    398     t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2)); 
    399  
    400     t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    401  
    402     t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    403  
    404     t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu)); 
    405  
    406     t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))); 
    407  
    408     t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2)))); 
    409  
    410     t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))); 
    411  
    412     t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu)))); 
    413  
    414  
    415     yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))))); 
    416  
    417     return (yy); 
    418 } 
    419  
    420  
    421 static double 
    422 Sk_WR(double q, double L, double b) 
    423 { 
    424         // 
    425         double p1,p2,p1short,p2short,q0; 
    426         double C,ans,q0short,Sexvmodify,pi; 
    427  
    428         pi = 4.0*atan(1.0); 
    429  
    430         p1 = 4.12; 
    431         p2 = 4.42; 
    432         p1short = 5.36; 
    433         p2short = 5.62; 
    434         q0 = 3.1; 
    435         // 
    436         q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0); 
    437  
    438         // 
    439         if(L/b > 10.0) { 
    440                 C = 3.06/pow((L/b),0.44); 
    441         } else { 
    442                 C = 1.0; 
    443         } 
    444         // 
    445  
    446         if( L > 4*b ) { // Longer Chains 
    447                 if (q*b <= 3.1) {               //Modified by Yun on Oct. 15, 
    448                         Sexvmodify = Sexvnew(q, L, b); 
    449                         ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L); 
    450                 } else { //q(i)*b > 3.1 
    451                         ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L); 
    452                 } 
    453         } else { //L <= 4*b Shorter Chains 
    454                 if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) { 
    455                         if (q*b<=0.01) { 
    456                                 ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0; 
    457                         } else { 
    458                                 ans = Sdebye1(q,L,b); 
    459                         } 
    460                 } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3) 
    461                         ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L); 
    462                 } 
    463         } 
    464  
    465         return(ans); 
    466         //return(a2long(q, L, b, p1, p2, q0)); 
    467 } 
    468  

  • sansmodels/src/sans/models/libigor/libSphere.c

    r34c2649 r5f0dcab  
    88#include "GaussWeights.h" 
    99#include "libSphere.h" 
     10 
     11 
     12static double 
     13gammln(double xx) { 
     14    double x,y,tmp,ser; 
     15    static double cof[6]={76.18009172947146,-86.50532032941677, 
     16    24.01409824083091,-1.231739572450155, 
     17    0.1208650973866179e-2,-0.5395239384953e-5}; 
     18    int j; 
     19 
     20    y=x=xx; 
     21    tmp=x+5.5; 
     22    tmp -= (x+0.5)*log(tmp); 
     23    ser=1.000000000190015; 
     24    for (j=0;j<=5;j++) ser += cof[j]/++y; 
     25    return -tmp+log(2.5066282746310005*ser/x); 
     26} 
     27 
     28static double 
     29LogNormal_distr(double sig, double mu, double pt) 
     30{ 
     31  double retval,pi; 
     32 
     33  pi = 4.0*atan(1.0); 
     34  retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
     35  return(retval); 
     36} 
     37 
     38static double 
     39Gauss_distr(double sig, double avg, double pt) 
     40{ 
     41  double retval,Pi; 
     42 
     43  Pi = 4.0*atan(1.0); 
     44  retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0); 
     45  return(retval); 
     46} 
     47 
     48static double SchulzPoint(double x, double avg, double zz) { 
     49    double dr; 
     50    dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); 
     51    return (exp(dr)); 
     52}; 
     53 
    1054 
    1155// scattering from a sphere - hardly needs to be an XOP... 

  • sansmodels/src/sans/models/libigor/libSphere.h

    r6a2c931 r5f0dcab  
    4949double SC_Integrand(double w[], double qq, double xx, double yy); 
    5050double SCeval(double Theta, double Phi, double temp3, double temp4, double temp5); 
    51  
    52  
    53  
    54 static double 
    55 gammln(double xx) { 
    56     double x,y,tmp,ser; 
    57     static double cof[6]={76.18009172947146,-86.50532032941677, 
    58                 24.01409824083091,-1.231739572450155, 
    59                 0.1208650973866179e-2,-0.5395239384953e-5}; 
    60     int j; 
    61  
    62     y=x=xx; 
    63     tmp=x+5.5; 
    64     tmp -= (x+0.5)*log(tmp); 
    65     ser=1.000000000190015; 
    66     for (j=0;j<=5;j++) ser += cof[j]/++y; 
    67     return -tmp+log(2.5066282746310005*ser/x); 
    68 } 
    69  
    70 static double 
    71 LogNormal_distr(double sig, double mu, double pt) 
    72 { 
    73         double retval,pi; 
    74  
    75         pi = 4.0*atan(1.0); 
    76         retval = (1.0/ (sig*pt*sqrt(2.0*pi)) )*exp( -0.5*(log(pt) - mu)*(log(pt) - mu)/sig/sig ); 
    77         return(retval); 
    78 } 
    79  
    80 static double 
    81 Gauss_distr(double sig, double avg, double pt) 
    82 { 
    83         double retval,Pi; 
    84  
    85         Pi = 4.0*atan(1.0); 
    86         retval = (1.0/ (sig*sqrt(2.0*Pi)) )*exp(-(avg-pt)*(avg-pt)/sig/sig/2.0); 
    87         return(retval); 
    88 } 
    89  
    90 static double SchulzPoint(double x, double avg, double zz) { 
    91     double dr; 
    92     dr = zz*log(x) - gammln(zz+1.0)+(zz+1.0)*log((zz+1.0)/avg)-(x/avg*(zz+1.0)); 
    93     return (exp(dr)); 
    94 }; 
    95  
Note: See TracChangeset for help on using the changeset viewer.