Changeset 34c2649 in sasview for sansmodels/src/sans/models/libigor/libCylinder.c

Timestamp:

Nov 21, 2011 2:52:25 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:

1810613

Parents:

afea8fe

Message:

Re #4 Fixed warnings

File:

• sansmodels/src/sans/models/libigor/libCylinder.c (modified) (2 diffs)

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

@@ -1337 +1337 @@
         
         return answer;
-}
-
-/////////functions for WRC implementation of flexible cylinders
-static double
-Sk_WR(double q, double L, double b)
-{
-        //
-        double p1,p2,p1short,p2short,q0,qconnect;
-        double C,epsilon,ans,q0short,Sexvmodify,pi;
-        
-        pi = 4.0*atan(1.0);
-        
-        p1 = 4.12;
-        p2 = 4.42;
-        p1short = 5.36;
-        p2short = 5.62;
-        q0 = 3.1;
-        qconnect = q0/b;
-        //      
-        q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);
-        
-        //
-        if(L/b > 10.0) {
-                C = 3.06/pow((L/b),0.44);
-                epsilon = 0.176;
-        } else {
-                C = 1.0;
-                epsilon = 0.170;
-        }
-        //
-        
-        if( L > 4*b ) { // Longer Chains
-                if (q*b <= 3.1) {               //Modified by Yun on Oct. 15,
-                        Sexvmodify = Sexvnew(q, L, b);
-                        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);
-                } else { //q(i)*b > 3.1
-                        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);
-                } 
-        } else { //L <= 4*b Shorter Chains
-                if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
-                        if (q*b<=0.01) {
-                                ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
-                        } else {
-                                ans = Sdebye1(q,L,b);
-                        }
-                } else {        //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
-                        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);
-                }
-        }
-        
-        return(ans);
-        //return(a2long(q, L, b, p1, p2, q0));
-}
-
-//WR named this w (too generic)
-static double
-w_WR(double x)
-{
-    double yy;
-    yy = 0.5*(1 + tanh((x - 1.523)/0.1477));
-        
-    return (yy);
-}
-
-//
-static double
-u1(double q, double L, double b)
-{
-    double yy;
-        
-    yy = Rgsquareshort(q,L,b)*q*q;
-    
-    return (yy);
-}
-
-// was named u
-static double
-u_WR(double q, double L, double b)
-{
-    double yy;
-    yy = Rgsquare(q,L,b)*q*q;
-    return (yy);
-}
-
-
-
-//
-static double
-Rgsquarezero(double q, double L, double b)
-{    
-    double yy;
-    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))));
-    
-    return (yy);
-}
-
-//
-static double
-Rgsquareshort(double q, double L, double b)
-{    
-    double yy;
-    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);
-    
-    return (yy);
-}
-
-//
-static double
-Rgsquare(double q, double L, double b)
-{
-    double yy;
-    yy = AlphaSquare(L/b)*L*b/6.0;
-    
-    return (yy);
-}
-
-//
-static double
-AlphaSquare(double x)
-{    
-    double yy;
-    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );
-        
-    return (yy);
-}
-
-// ?? funciton is not used - but should the log actually be log10???
-static double
-miu(double x)
-{    
-    double yy;
-    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)));
-    
-    return (yy);
-}
-
-//
-static double
-Sdebye(double q, double L, double b)
-{    
-    double yy;
-    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));
-        
-    return (yy);
-}
-
-//
-static double
-Sdebye1(double q, double L, double b)
-{    
-    double yy;
-    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );
-    
-    return (yy);
-}
-
-//
-static double
-Sexv(double q, double L, double b)
-{    
-    double yy,C1,C2,C3,miu,Rg2;
-    C1=1.22;
-    C2=0.4288;
-    C3=-1.651;
-    miu = 0.585;
-        
-    Rg2 = Rgsquare(q,L,b);
-    
-    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)));
-    
-    return (yy);
-}
-
-// this must be WR modified version
-static double
-Sexvnew(double q, double L, double b)
-{    
-    double yy,C1,C2,C3,miu;
-    double del=1.05,C_star2,Rg2;
-    
-    C1=1.22;
-    C2=0.4288;
-    C3=-1.651;
-    miu = 0.585;
-        
-    //calculating the derivative to decide on the corection (cutoff) term?
-    // I have modified this from WRs original code
-    
-    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
-        C_star2 = 0.0;
-    } else {
-        C_star2 = 1.0;
-    }
-        
-    Rg2 = Rgsquare(q,L,b);
-    
-        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)));
-        
-    return (yy);
-}
-
-// these are the messy ones
-static double
-a2short(double q, double L, double b, double p1short, double p2short, double q0)
-{    
-    double yy,Rg2_sh;
-    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
-    double pi;
-        
-    E = 2.718281828459045091;
-        pi = 4.0*atan(1.0);
-    Rg2_sh = Rgsquareshort(q,L,b);
-    Rg2_sh2 = Rg2_sh*Rg2_sh;
-    t1 = ((q0*q0*Rg2_sh)/(b*b));
-    Et1 = pow(E,t1);
-    Emt1 =pow(E,-t1); 
-    q02 = q0*q0;
-    q0p = pow(q0,(-4.0 + p2short) );
-    
-    //E is the number e
-    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)))))));
-        
-    return (yy);
-}
-
-//
-static double
-a1short(double q, double L, double b, double p1short, double p2short, double q0)
-{    
-    double yy,Rg2_sh;
-    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
-        double pi;
-    
-    E = 2.718281828459045091;
-        pi = 4.0*atan(1.0);
-    Rg2_sh = Rgsquareshort(q,L,b);
-    Rg2_sh2 = Rg2_sh*Rg2_sh;
-    b3 = b*b*b;
-    t1 = ((q0*q0*Rg2_sh)/(b*b));
-    Et1 = pow(E,t1);
-    Emt1 =pow(E,-t1); 
-    q02 = q0*q0;
-    q0p = pow(q0,(-4.0 + p1short) );
-
-    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))))));
-
-    return(yy);
-}
-
-// this one will be lots of trouble
-static double
-a2long(double q, double L, double b, double p1, double p2, double q0)
-{
-    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
-    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
-    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;
-    
-    pi = 4.0*atan(1.0);
-    E = 2.718281828459045091;
-    if( L/b > 10.0) {
-        C = 3.06/pow((L/b),0.44);
-    } else {
-        C = 1.0;
-    }
-        
-    C1 = 1.22;
-    C2 = 0.4288;
-    C3 = -1.651;
-    C4 = 1.523;
-    C5 = 0.1477;
-    miu = 0.585;
-        
-    Rg2 = Rgsquare(q,L,b);
-    Rg22 = Rg2*Rg2;
-    b2 = b*b;
-    b3 = b*b*b;
-    b4 = b3*b;
-    q02 = q0*q0;
-    q03 = q0*q0*q0;
-    q04 = q03*q0;
-    q05 = q04*q0;
-    
-    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));
-
-    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;
-
-    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);
-
-    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);
-
-    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);
-
-    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);
-
-    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));
-
-    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-
-    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;
-
-    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);
-
-
-    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))))))))));
-
-    return (yy);
-}
-
-//need to define this on my own
-static double
-sech_WR(double x)
-{       
-        return(1/cosh(x));
-}
-
-//
-static double
-a1long(double q, double L, double b, double p1, double p2, double q0)
-{    
-    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
-    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
-    double E,pi;
-    double b2,b3,b4,q02,q03,q04,q05,Rg22;
-        
-    pi = 4.0*atan(1.0);
-    E = 2.718281828459045091;
-    
-    if( L/b > 10.0) {
-        C = 3.06/pow((L/b),0.44);
-    } else {
-        C = 1.0;
-    }
-        
-    C1 = 1.22;
-    C2 = 0.4288;
-    C3 = -1.651;
-    C4 = 1.523;
-    C5 = 0.1477;
-    miu = 0.585;
-        
-    Rg2 = Rgsquare(q,L,b);
-    Rg22 = Rg2*Rg2;
-    b2 = b*b;
-    b3 = b*b*b;
-    b4 = b3*b;
-    q02 = q0*q0;
-    q03 = q0*q0*q0;
-    q04 = q03*q0;
-    q05 = q04*q0;
-    
-    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))));
-    
-    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))))));
-    
-    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))));
-    
-    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-        
-    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));
-    
-    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)));
-    
-    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));
-    
-    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));
-    
-    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))))));
-    
-    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))))));
-    
-    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));
-    
-    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));
-    
-    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))));
-    
-    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))))));
-    
-    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))));
-        
-    
-    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)))))))))));
-        
-    return (yy);
 }
 
@@ -2132 +1748 @@
 }
 
-static double
-gammaln(double xx)
-{
-        double x,y,tmp,ser;
-        static double cof[6]={76.18009172947146,-86.50532032941677,
-                24.01409824083091,-1.231739572450155,
-                0.1208650973866179e-2,-0.5395239384953e-5};
-        int j;
-        
-        y=x=xx;
-        tmp=x+5.5;
-        tmp -= (x+0.5)*log(tmp);
-        ser=1.000000000190015;
-        for (j=0;j<=5;j++) ser += cof[j]/++y;
-        return -tmp+log(2.5066282746310005*ser/x);
-}
-
 
 double