Changeset a629d8e in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 11:40:06 AM (9 years ago)
Author:
wojciech
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
3a45c2c, 95ce773
Parents:
e17904a (diff), 292fe08 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of ​https://github.com/SasView/sasmodels

Location:
sasmodels/models
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/hardsphere.py

    re98c1e0 r934f906  
    7575Iq = """ 
    7676      double D,A,B,G,X,X2,X4,S,C,FF,HARDSPH; 
     77      // these are c compiler instructions, can also put normal code inside the "if else" structure  
     78      #if FLOAT_SIZE > 4 
     79      // double precision    orig had 0.2, don't call the variable cutoff as PAK already has one called that! Must use UPPERCASE name please. 
     80      //  0.05 better, 0.1 OK 
     81      #define CUTOFFHS 0.05 
     82      #else 
     83      // 0.1 bad, 0.2 OK, 0.3 good, 0.4 better, 0.8 no good 
     84      #define CUTOFFHS 0.4   
     85      #endif 
    7786 
    7887      if(fabs(radius_effective) < 1.E-12) { 
    … …  
    97106      G=0.5*volfraction*A; 
    98107 
    99       if(X < 0.2) { 
    100       // RKH Feb 2016, use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double.  
     108      if(X < CUTOFFHS) { 
     109      // RKH Feb 2016, use Taylor series expansion for small X 
    101110      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
    102111      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient  

  • sasmodels/models/lib/j1_cephes.c

    ra5af4e1 rfad5dc1  
    3939Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    4040*/ 
    41 double J1(double ); 
     41 
     42constant double RP[8] = { 
     43    -8.99971225705559398224E8, 
     44    4.52228297998194034323E11, 
     45    -7.27494245221818276015E13, 
     46    3.68295732863852883286E15, 
     47    0.0, 
     48    0.0, 
     49    0.0, 
     50    0.0 }; 
     51 
     52constant double RQ[8] = { 
     53        6.20836478118054335476E2, 
     54        2.56987256757748830383E5, 
     55        8.35146791431949253037E7, 
     56        2.21511595479792499675E10, 
     57        4.74914122079991414898E12, 
     58        7.84369607876235854894E14, 
     59        8.95222336184627338078E16, 
     60        5.32278620332680085395E18 
     61    }; 
     62 
     63constant double PP[8] = { 
     64    7.62125616208173112003E-4, 
     65    7.31397056940917570436E-2, 
     66    1.12719608129684925192E0, 
     67    5.11207951146807644818E0, 
     68    8.42404590141772420927E0, 
     69    5.21451598682361504063E0, 
     70    1.00000000000000000254E0, 
     71    0.0} ; 
     72 
     73 
     74constant double PQ[8] = { 
     75    5.71323128072548699714E-4, 
     76    6.88455908754495404082E-2, 
     77    1.10514232634061696926E0, 
     78    5.07386386128601488557E0, 
     79    8.39985554327604159757E0, 
     80    5.20982848682361821619E0, 
     81    9.99999999999999997461E-1, 
     82    0.0 }; 
     83 
     84constant double QP[8] = { 
     85    5.10862594750176621635E-2, 
     86    4.98213872951233449420E0, 
     87    7.58238284132545283818E1, 
     88    3.66779609360150777800E2, 
     89    7.10856304998926107277E2, 
     90    5.97489612400613639965E2, 
     91    2.11688757100572135698E2, 
     92    2.52070205858023719784E1 }; 
     93 
     94constant double QQ[8] = { 
     95    7.42373277035675149943E1, 
     96    1.05644886038262816351E3, 
     97    4.98641058337653607651E3, 
     98    9.56231892404756170795E3, 
     99    7.99704160447350683650E3, 
     100    2.82619278517639096600E3, 
     101    3.36093607810698293419E2, 
     102    0.0 }; 
     103 
     104constant double JP[8] = { 
     105        -4.878788132172128E-009, 
     106        6.009061827883699E-007, 
     107        -4.541343896997497E-005, 
     108        1.937383947804541E-003, 
     109        -3.405537384615824E-002, 
     110        0.0, 
     111        0.0, 
     112        0.0 
     113    }; 
     114 
     115constant double MO1[8] = { 
     116        6.913942741265801E-002, 
     117        -2.284801500053359E-001, 
     118        3.138238455499697E-001, 
     119        -2.102302420403875E-001, 
     120        5.435364690523026E-003, 
     121        1.493389585089498E-001, 
     122        4.976029650847191E-006, 
     123        7.978845453073848E-001 
     124    }; 
     125 
     126constant double PH1[8] = { 
     127        -4.497014141919556E+001, 
     128        5.073465654089319E+001, 
     129        -2.485774108720340E+001, 
     130        7.222973196770240E+000, 
     131        -1.544842782180211E+000, 
     132        3.503787691653334E-001, 
     133        -1.637986776941202E-001, 
     134        3.749989509080821E-001 
     135    }; 
    42136 
    43137double J1(double x) { 
    … …  
    60154        { 
    61155            z = x * x; 
    62             w = polevlRP( z, 3 ) / p1evlRQ( z, 8 ); 
     156            w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 
    63157            w = w * x * (z - Z1) * (z - Z2); 
    64158            return( w ); 
    … …  
    68162    z = w * w; 
    69163 
    70     p = polevlPP( z, 6)/polevlPQ( z, 6 ); 
    71     q = polevlQP( z, 7)/p1evlQQ( z, 7 ); 
     164    p = polevl( z, PP, 6)/polevl( z, PQ, 6 ); 
     165    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 ); 
    72166 
    73167    xn = x - THPIO4; 
    … …  
    95189        { 
    96190            z = xx * xx; 
    97             p = (z-Z1) * xx * polevlJP( z, 4 ); 
     191            p = (z-Z1) * xx * polevl( z, JP, 4 ); 
    98192            return( p ); 
    99193        } 
    … …  
    102196    w = sqrt(q); 
    103197 
    104     p = w * polevlMO1( q, 7); 
     198    p = w * polevl( q, MO1, 7); 
    105199    w = q*q; 
    106     xn = q * polevlPH1( w, 7) - THPIO4F; 
     200    xn = q * polevl( w, PH1, 7) - THPIO4F; 
    107201    p = p * cos(xn + xx); 
    108202 
    … …  
    110204#endif 
    111205} 
     206 
     207 
     208/*double J1c(double x) { 
     209    return 
     210}*/ 

  • sasmodels/models/lib/polevl.c

    r3b12dea rfad5dc1  
    5151*/ 
    5252 
    53 double polevlRP(double x, int N ) { 
    5453 
    55     double coef[8] = { 
    56     -8.99971225705559398224E8, 
    57     4.52228297998194034323E11, 
    58     -7.27494245221818276015E13, 
    59     3.68295732863852883286E15, 
    60     0.0, 
    61     0.0, 
    62     0.0, 
    63     0.0 }; 
    64  
    65     int i = 0; 
    66     double ans = coef[i]; 
    67  
    68     while (i < N) { 
    69         i++; 
    70         ans = ans * x + coef[i]; 
    71     } 
    72     return ans ; 
    73 } 
    74  
    75 double polevlRQ(double x, int N ) { 
    76  
    77     double coef[8] = { 
    78         6.20836478118054335476E2, 
    79         2.56987256757748830383E5, 
    80         8.35146791431949253037E7, 
    81         2.21511595479792499675E10, 
    82         4.74914122079991414898E12, 
    83         7.84369607876235854894E14, 
    84         8.95222336184627338078E16, 
    85         5.32278620332680085395E18 
    86     }; 
    87  
    88     int i = 0; 
    89     double ans = coef[i]; 
    90  
    91     while (i < N) { 
    92         i++; 
    93         ans = ans * x + coef[i]; 
    94     } 
    95     return ans ; 
    96 } 
    97  
    98 double polevlPP(double x, int N ) { 
    99  
    100     double coef[8] = { 
    101     7.62125616208173112003E-4, 
    102     7.31397056940917570436E-2, 
    103     1.12719608129684925192E0, 
    104     5.11207951146807644818E0, 
    105     8.42404590141772420927E0, 
    106     5.21451598682361504063E0, 
    107     1.00000000000000000254E0, 
    108     0.0} ; 
    109  
    110     int i = 0; 
    111     double ans = coef[i]; 
    112  
    113     while (i < N) { 
    114         i++; 
    115         ans = ans * x + coef[i]; 
    116     } 
    117     return ans ; 
    118 } 
    119  
    120 double polevlPQ(double x, int N ) { 
    121  
    122     double coef[8] = { 
    123     5.71323128072548699714E-4, 
    124     6.88455908754495404082E-2, 
    125     1.10514232634061696926E0, 
    126     5.07386386128601488557E0, 
    127     8.39985554327604159757E0, 
    128     5.20982848682361821619E0, 
    129     9.99999999999999997461E-1, 
    130     0.0 }; 
    131  
    132     int i = 0; 
    133     double ans = coef[i]; 
    134  
    135     while (i < N) { 
    136         i++; 
    137         ans = ans * x + coef[i]; 
    138     } 
    139     return ans ; 
    140 } 
    141  
    142 double polevlQP(double x, int N ) { 
    143  
    144     double coef[8] = { 
    145     5.10862594750176621635E-2, 
    146     4.98213872951233449420E0, 
    147     7.58238284132545283818E1, 
    148     3.66779609360150777800E2, 
    149     7.10856304998926107277E2, 
    150     5.97489612400613639965E2, 
    151     2.11688757100572135698E2, 
    152     2.52070205858023719784E1 }; 
    153  
    154     int i = 0; 
    155     double ans = coef[i]; 
    156  
    157     while (i < N) { 
    158         i++; 
    159         ans = ans * x + coef[i]; 
    160     } 
    161     return ans ; 
    162  
    163 } 
    164  
    165 double polevlQQ(double x, int N ) { 
    166  
    167     double coef[8] = { 
    168     7.42373277035675149943E1, 
    169     1.05644886038262816351E3, 
    170     4.98641058337653607651E3, 
    171     9.56231892404756170795E3, 
    172     7.99704160447350683650E3, 
    173     2.82619278517639096600E3, 
    174     3.36093607810698293419E2, 
    175     0.0 }; 
    176  
    177     int i = 0; 
    178     double ans = coef[i]; 
    179  
    180     while (i < N) { 
    181         i++; 
    182         ans = ans * x + coef[i]; 
    183     } 
    184     return ans ; 
    185  
    186 } 
    187  
    188 double polevlJP( double x, int N ) { 
    189     double coef[8] = { 
    190         -4.878788132172128E-009, 
    191         6.009061827883699E-007, 
    192         -4.541343896997497E-005, 
    193         1.937383947804541E-003, 
    194         -3.405537384615824E-002, 
    195         0.0, 
    196         0.0, 
    197         0.0 
    198     }; 
    199  
    200     int i = 0; 
    201     double ans = coef[i]; 
    202  
    203     while (i < N) { 
    204         i++; 
    205         ans = ans * x + coef[i]; 
    206     } 
    207     return ans ; 
    208  
    209 } 
    210  
    211 double polevlMO1( double x, int N ) { 
    212     double coef[8] = { 
    213         6.913942741265801E-002, 
    214         -2.284801500053359E-001, 
    215         3.138238455499697E-001, 
    216         -2.102302420403875E-001, 
    217         5.435364690523026E-003, 
    218         1.493389585089498E-001, 
    219         4.976029650847191E-006, 
    220         7.978845453073848E-001 
    221     }; 
    222  
    223     int i = 0; 
    224     double ans = coef[i]; 
    225  
    226     while (i < N) { 
    227         i++; 
    228         ans = ans * x + coef[i]; 
    229     } 
    230     return ans ; 
    231 } 
    232  
    233 double polevlPH1( double x, int N ) { 
    234  
    235     double coef[8] = { 
    236         -4.497014141919556E+001, 
    237         5.073465654089319E+001, 
    238         -2.485774108720340E+001, 
    239         7.222973196770240E+000, 
    240         -1.544842782180211E+000, 
    241         3.503787691653334E-001, 
    242         -1.637986776941202E-001, 
    243         3.749989509080821E-001 
    244     }; 
    245  
    246     int i = 0; 
    247     double ans = coef[i]; 
    248  
    249     while (i < N) { 
    250         i++; 
    251         ans = ans * x + coef[i]; 
    252     } 
    253     return ans ; 
    254 } 
    255  
    256 /*double polevl( double x, double coef[8], int N ) { 
     54double polevl( double x, constant double *coef, int N ) { 
    25755 
    25856    int i = 0; 
    … …  
    26664    return ans ; 
    26765 
    268 }*/ 
     66} 
    26967 
    27068/*                                                      p1evl() */ 
    … …  
    27472 */ 
    27573 
    276 double p1evlRP( double x, int N ) { 
    277     double coef[8] = { 
    278     -8.99971225705559398224E8, 
    279     4.52228297998194034323E11, 
    280     -7.27494245221818276015E13, 
    281     3.68295732863852883286E15, 
    282     0.0, 
    283     0.0, 
    284     0.0, 
    285     0.0 }; 
    286  
    287     int i=0; 
    288     double ans = x+coef[i]; 
    289  
    290     while (i < N-1) { 
    291         i++; 
    292         ans = ans*x + coef[i]; 
    293     } 
    294  
    295     return( ans ); 
    296  
    297 } 
    298  
    299 double p1evlRQ( double x, int N ) { 
    300 //1: RQ 
    301     double coef[8] = { 
    302     6.20836478118054335476E2, 
    303     2.56987256757748830383E5, 
    304     8.35146791431949253037E7, 
    305     2.21511595479792499675E10, 
    306     4.74914122079991414898E12, 
    307     7.84369607876235854894E14, 
    308     8.95222336184627338078E16, 
    309     5.32278620332680085395E18}; 
    310  
     74double p1evl( double x, constant double *coef, int N ) { 
    31175    int i=0; 
    31276    double ans = x+coef[i]; 
    … …  
    31983    return( ans ); 
    32084} 
    321  
    322 double p1evlPP( double x, int N ) { 
    323 //3 : PP 
    324     double coef[8] = { 
    325     7.62125616208173112003E-4, 
    326     7.31397056940917570436E-2, 
    327     1.12719608129684925192E0, 
    328     5.11207951146807644818E0, 
    329     8.42404590141772420927E0, 
    330     5.21451598682361504063E0, 
    331     1.00000000000000000254E0, 
    332     0.0}; 
    333  
    334     int i=0; 
    335     double ans = x+coef[i]; 
    336  
    337     while (i < N-1) { 
    338         i++; 
    339         ans = ans*x + coef[i]; 
    340     } 
    341  
    342     return( ans ); 
    343 } 
    344  
    345 double p1evlPQ( double x, int N ) { 
    346 //4: PQ 
    347     double coef[8] = { 
    348     5.71323128072548699714E-4, 
    349     6.88455908754495404082E-2, 
    350     1.10514232634061696926E0, 
    351     5.07386386128601488557E0, 
    352     8.39985554327604159757E0, 
    353     5.20982848682361821619E0, 
    354     9.99999999999999997461E-1, 
    355     0.0}; 
    356  
    357     int i=0; 
    358     double ans = x+coef[i]; 
    359  
    360     while (i < N-1) { 
    361         i++; 
    362         ans = ans*x + coef[i]; 
    363     } 
    364  
    365     return( ans ); 
    366 } 
    367  
    368 double p1evlQP( double x, int N ) { 
    369 //5: QP 
    370     double coef[8] = { 
    371     5.10862594750176621635E-2, 
    372     4.98213872951233449420E0, 
    373     7.58238284132545283818E1, 
    374     3.66779609360150777800E2, 
    375     7.10856304998926107277E2, 
    376     5.97489612400613639965E2, 
    377     2.11688757100572135698E2, 
    378     2.52070205858023719784E1 }; 
    379  
    380     int i=0; 
    381     double ans = x+coef[i]; 
    382  
    383     while (i < N-1) { 
    384         i++; 
    385         ans = ans*x + coef[i]; 
    386     } 
    387  
    388     return( ans ); 
    389 } 
    390  
    391 double p1evlQQ( double x, int N ) { 
    392     double coef[8] = { 
    393     7.42373277035675149943E1, 
    394     1.05644886038262816351E3, 
    395     4.98641058337653607651E3, 
    396     9.56231892404756170795E3, 
    397     7.99704160447350683650E3, 
    398     2.82619278517639096600E3, 
    399     3.36093607810698293419E2, 
    400     0.0}; 
    401  
    402     int i=0; 
    403     double ans = x+coef[i]; 
    404  
    405     while (i < N-1) { 
    406         i++; 
    407         ans = ans*x + coef[i]; 
    408     } 
    409  
    410     return( ans ); 
    411  
    412 } 
    413  
    414 double p1evlJP( double x, int N ) { 
    415     double coef[8] = { 
    416         -4.878788132172128E-009, 
    417         6.009061827883699E-007, 
    418         -4.541343896997497E-005, 
    419         1.937383947804541E-003, 
    420         -3.405537384615824E-002, 
    421         0.0, 
    422         0.0, 
    423         0.0}; 
    424  
    425     int i=0; 
    426     double ans = x+coef[i]; 
    427  
    428     while (i < N-1) { 
    429         i++; 
    430         ans = ans*x + coef[i]; 
    431     } 
    432  
    433     return( ans ); 
    434 } 
    435  
    436 double p1evlMO1( double x, int N ) { 
    437     double coef[8] = { 
    438         6.913942741265801E-002, 
    439         -2.284801500053359E-001, 
    440         3.138238455499697E-001, 
    441         -2.102302420403875E-001, 
    442         5.435364690523026E-003, 
    443         1.493389585089498E-001, 
    444         4.976029650847191E-006, 
    445         7.978845453073848E-001}; 
    446  
    447     int i=0; 
    448     double ans = x+coef[i]; 
    449  
    450     while (i < N-1) { 
    451         i++; 
    452         ans = ans*x + coef[i]; 
    453     } 
    454  
    455     return( ans ); 
    456  
    457 } 
    458  
    459 double p1evlPH1( double x, int N ) { 
    460     double coef[8] = { 
    461         -4.497014141919556E+001, 
    462         5.073465654089319E+001, 
    463         -2.485774108720340E+001, 
    464         7.222973196770240E+000, 
    465         -1.544842782180211E+000, 
    466         3.503787691653334E-001, 
    467         -1.637986776941202E-001, 
    468         3.749989509080821E-001}; 
    469     int i=0; 
    470     double ans = x+coef[i]; 
    471  
    472     while (i < N-1) { 
    473         i++; 
    474         ans = ans*x + coef[i]; 
    475     } 
    476  
    477     return( ans ); 
    478 } 
    479  
    480 /*double p1evl( double x, double coef[8], int N ) { 
    481     int i=0; 
    482     double ans = x+coef[i]; 
    483  
    484     while (i < N-1) { 
    485         i++; 
    486         ans = ans*x + coef[i]; 
    487     } 
    488  
    489     return( ans ); 
    490  
    491 }*/ 
Note: See TracChangeset for help on using the changeset viewer.