Changeset a5af4e1 in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 10:39:55 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:
7904790
Parents:
c094758
Message:

Polevl fixed

Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/bessel.py

    rcbd37a7 ra5af4e1  
    7878 
    7979Iq = """ 
    80     return 2.0*j1(q)/q; 
     80    return J1(q); 
    8181    """ 
    8282 

  • sasmodels/models/lib/j1_cephes.c

    re2af2a9 ra5af4e1  
    3939Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    4040*/ 
    41 double j1(double ); 
     41double J1(double ); 
    4242 
    43 double j1(double x) { 
     43double J1(double x) { 
    4444 
    4545//Cephes double pression function 
    … …  
    4848    double w, z, p, q, xn; 
    4949 
    50     const double DR1 = 5.78318596294678452118E0; 
    51     const double DR2 = 3.04712623436620863991E1; 
    5250    const double Z1 = 1.46819706421238932572E1; 
    5351    const double Z2 = 4.92184563216946036703E1; 
    5452    const double THPIO4 =  2.35619449019234492885; 
    5553    const double SQ2OPI = 0.79788456080286535588; 
    56  
    57     double RP[8] = { 
    58     -8.99971225705559398224E8, 
    59     4.52228297998194034323E11, 
    60     -7.27494245221818276015E13, 
    61     3.68295732863852883286E15, 
    62     0.0, 
    63     0.0, 
    64     0.0, 
    65     0.0 
    66     }; 
    67  
    68     double RQ[8] = { 
    69     /* 1.00000000000000000000E0,*/ 
    70     6.20836478118054335476E2, 
    71     2.56987256757748830383E5, 
    72     8.35146791431949253037E7, 
    73     2.21511595479792499675E10, 
    74     4.74914122079991414898E12, 
    75     7.84369607876235854894E14, 
    76     8.95222336184627338078E16, 
    77     5.32278620332680085395E18, 
    78     }; 
    79  
    80     double PP[8] = { 
    81     7.62125616208173112003E-4, 
    82     7.31397056940917570436E-2, 
    83     1.12719608129684925192E0, 
    84     5.11207951146807644818E0, 
    85     8.42404590141772420927E0, 
    86     5.21451598682361504063E0, 
    87     1.00000000000000000254E0, 
    88     0.0, 
    89     }; 
    90     double PQ[8] = { 
    91     5.71323128072548699714E-4, 
    92     6.88455908754495404082E-2, 
    93     1.10514232634061696926E0, 
    94     5.07386386128601488557E0, 
    95     8.39985554327604159757E0, 
    96     5.20982848682361821619E0, 
    97     9.99999999999999997461E-1, 
    98     0.0, 
    99     }; 
    100  
    101     double QP[8] = { 
    102     5.10862594750176621635E-2, 
    103     4.98213872951233449420E0, 
    104     7.58238284132545283818E1, 
    105     3.66779609360150777800E2, 
    106     7.10856304998926107277E2, 
    107     5.97489612400613639965E2, 
    108     2.11688757100572135698E2, 
    109     2.52070205858023719784E1, 
    110     }; 
    111  
    112     double QQ[8] = { 
    113     /* 1.00000000000000000000E0,*/ 
    114     7.42373277035675149943E1, 
    115     1.05644886038262816351E3, 
    116     4.98641058337653607651E3, 
    117     9.56231892404756170795E3, 
    118     7.99704160447350683650E3, 
    119     2.82619278517639096600E3, 
    120     3.36093607810698293419E2, 
    121     0.0, 
    122     }; 
    12354 
    12455    w = x; 
    … …  
    12960        { 
    13061            z = x * x; 
    131             w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 
     62            w = polevlRP( z, 3 ) / p1evlRQ( z, 8 ); 
    13263            w = w * x * (z - Z1) * (z - Z2); 
    13364            return( w ); 
    … …  
    13667    w = 5.0/x; 
    13768    z = w * w; 
    138     p = polevl( z, PP, 6)/polevl( z, PQ, 6 ); 
    139     q = polevl( z, QP, 7)/p1evl( z, QQ, 7 ); 
     69 
     70    p = polevlPP( z, 6)/polevlPQ( z, 6 ); 
     71    q = polevlQP( z, 7)/p1evlQQ( z, 7 ); 
     72 
    14073    xn = x - THPIO4; 
    14174 
    … …  
    15588 
    15689 
    157     double JP[8] = { 
    158         -4.878788132172128E-009, 
    159         6.009061827883699E-007, 
    160         -4.541343896997497E-005, 
    161         1.937383947804541E-003, 
    162         -3.405537384615824E-002, 
    163         0.0, 
    164         0.0, 
    165         0.0 
    166     }; 
    167  
    168     double MO1[8] = { 
    169         6.913942741265801E-002, 
    170         -2.284801500053359E-001, 
    171         3.138238455499697E-001, 
    172         -2.102302420403875E-001, 
    173         5.435364690523026E-003, 
    174         1.493389585089498E-001, 
    175         4.976029650847191E-006, 
    176         7.978845453073848E-001 
    177     }; 
    178  
    179     double PH1[8] = { 
    180         -4.497014141919556E+001, 
    181         5.073465654089319E+001, 
    182         -2.485774108720340E+001, 
    183         7.222973196770240E+000, 
    184         -1.544842782180211E+000, 
    185         3.503787691653334E-001, 
    186         -1.637986776941202E-001, 
    187         3.749989509080821E-001 
    188     }; 
    189  
    19090    xx = x; 
    19191    if( xx < 0 ) 
    … …  
    19595        { 
    19696            z = xx * xx; 
    197             p = (z-Z1) * xx * polevl( z, JP, 4 ); 
     97            p = (z-Z1) * xx * polevlJP( z, 4 ); 
    19898            return( p ); 
    19999        } 
    … …  
    202102    w = sqrt(q); 
    203103 
    204     p = w * polevl( q, MO1, 7); 
     104    p = w * polevlMO1( q, 7); 
    205105    w = q*q; 
    206     xn = q * polevl( w, PH1, 7) - THPIO4F; 
     106    xn = q * polevlPH1( w, 7) - THPIO4F; 
    207107    p = p * cos(xn + xx); 
    208108 
    … …  
    210110#endif 
    211111} 
    212  

  • sasmodels/models/lib/polevl.c

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