Changes in / [68b8734:de97440] in sasmodels


Ignore:
Files:
7 added
10 deleted
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/bessel.py

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

  • sasmodels/models/lib/j1_cephes.c

    ra5af4e1 re2af2a9  
    3939Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier 
    4040*/ 
    41 double J1(double ); 
    42  
    43 double J1(double x) { 
     41double j1(double ); 
     42 
     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; 
    5052    const double Z1 = 1.46819706421238932572E1; 
    5153    const double Z2 = 4.92184563216946036703E1; 
    … …  
    5355    const double SQ2OPI = 0.79788456080286535588; 
    5456 
     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    }; 
     123 
    55124    w = x; 
    56125    if( x < 0 ) 
    … …  
    60129        { 
    61130            z = x * x; 
    62             w = polevlRP( z, 3 ) / p1evlRQ( z, 8 ); 
     131            w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 
    63132            w = w * x * (z - Z1) * (z - Z2); 
    64133            return( w ); 
    … …  
    67136    w = 5.0/x; 
    68137    z = w * w; 
    69  
    70     p = polevlPP( z, 6)/polevlPQ( z, 6 ); 
    71     q = polevlQP( z, 7)/p1evlQQ( z, 7 ); 
    72  
     138    p = polevl( z, PP, 6)/polevl( z, PQ, 6 ); 
     139    q = polevl( z, QP, 7)/p1evl( z, QQ, 7 ); 
    73140    xn = x - THPIO4; 
    74141 
    … …  
    88155 
    89156 
     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 
    90190    xx = x; 
    91191    if( xx < 0 ) 
    … …  
    95195        { 
    96196            z = xx * xx; 
    97             p = (z-Z1) * xx * polevlJP( z, 4 ); 
     197            p = (z-Z1) * xx * polevl( z, JP, 4 ); 
    98198            return( p ); 
    99199        } 
    … …  
    102202    w = sqrt(q); 
    103203 
    104     p = w * polevlMO1( q, 7); 
     204    p = w * polevl( q, MO1, 7); 
    105205    w = q*q; 
    106     xn = q * polevlPH1( w, 7) - THPIO4F; 
     206    xn = q * polevl( w, PH1, 7) - THPIO4F; 
    107207    p = p * cos(xn + xx); 
    108208 
    … …  
    110210#endif 
    111211} 
     212 

  • sasmodels/models/lib/polevl.c

    r3b12dea re2af2a9  
    5151*/ 
    5252 
    53 double polevlRP(double x, int N ) { 
     53double polevl( double x, double coef[8], int N ); 
     54double p1evl( double x, double coef[8], int N ); 
    5455 
    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 ) { 
     56double polevl( double x, double coef[8], int N ) { 
    25757 
    25858    int i = 0; 
    … …  
    26666    return ans ; 
    26767 
    268 }*/ 
     68} 
    26969 
    27070/*                                                      p1evl() */ 
    … …  
    27474 */ 
    27575 
    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  
     76double p1evl( double x, double coef[8], int N ) { 
    28777    int i=0; 
    28878    double ans = x+coef[i]; 
    … …  
    29686 
    29787} 
    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  
    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  
    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.