Changes in / [de97440:68b8734] in sasmodels


Ignore:
Files:
10 added
7 deleted
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 r3b12dea  
    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 
     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 
     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 
     165double 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 
     188double 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 
     211double 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 
     233double 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 ) { 
     257 
     258    int i = 0; 
     259    double ans = coef[i]; 
     260 
     261    while (i < N) { 
     262        i++; 
     263        ans = ans * x + coef[i]; 
     264    } 
     265 
     266    return ans ; 
     267 
     268}*/ 
    69269 
    70270/*                                                      p1evl() */ 
    … …  
    74274 */ 
    75275 
    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 } 
     276double 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 
     299double 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 
     311    int i=0; 
     312    double ans = x+coef[i]; 
     313 
     314    while (i < N-1) { 
     315        i++; 
     316        ans = ans*x + coef[i]; 
     317    } 
     318 
     319    return( ans ); 
     320} 
     321 
     322double 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 
     345double 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 
     368double 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 
     391double 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 
     414double 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 
     436double 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 
     459double 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.