00001 #include "cdflib.h"
00002 double erfc1(int *ind,double *x)
00003 /*
00004 -----------------------------------------------------------------------
00005          EVALUATION OF THE COMPLEMENTARY ERROR FUNCTION
00006  
00007           ERFC1(IND,X) = ERFC(X)            IF IND = 0
00008           ERFC1(IND,X) = EXP(X*X)*ERFC(X)   OTHERWISE
00009 -----------------------------------------------------------------------
00010 */
00011 {
00012 static double c = .564189583547756e0;
00013 static double a[5] = {
00014     .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
00015     .479137145607681e-01,.128379167095513e+00
00016 };
00017 static double b[3] = {
00018     .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
00019 };
00020 static double p[8] = {
00021     -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
00022     4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
00023     4.51918953711873e+02,3.00459261020162e+02
00024 };
00025 static double q[8] = {
00026     1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
00027     2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
00028     7.90950925327898e+02,3.00459260956983e+02
00029 };
00030 static double r[5] = {
00031     2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
00032     4.65807828718470e+00,2.82094791773523e-01
00033 };
00034 static double s[4] = {
00035     9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
00036     1.80124575948747e+01
00037 };
00038 static int K1 = 1;
00039 static double erfc1,ax,bot,e,t,top,w;
00040 /*
00041      ..
00042      .. Executable Statements ..
00043 */
00044 /*
00045                      ABS(X) .LE. 0.5
00046 */
00047     ax = fabs(*x);
00048     if(ax > 0.5e0) goto S10;
00049     t = *x**x;
00050     top = (((a[0]*t+a[1])*t+a[2])*t+a[3])*t+a[4]+1.0e0;
00051     bot = ((b[0]*t+b[1])*t+b[2])*t+1.0e0;
00052     erfc1 = 0.5e0+(0.5e0-*x*(top/bot));
00053     if(*ind != 0) erfc1 = exp(t)*erfc1;
00054     return erfc1;
00055 S10:
00056 /*
00057                   0.5 .LT. ABS(X) .LE. 4
00058 */
00059     if(ax > 4.0e0) goto S20;
00060     top = ((((((p[0]*ax+p[1])*ax+p[2])*ax+p[3])*ax+p[4])*ax+p[5])*ax+p[6])*ax+p[
00061       7];
00062     bot = ((((((q[0]*ax+q[1])*ax+q[2])*ax+q[3])*ax+q[4])*ax+q[5])*ax+q[6])*ax+q[
00063       7];
00064     erfc1 = top/bot;
00065     goto S40;
00066 S20:
00067 /*
00068                       ABS(X) .GT. 4
00069 */
00070     if(*x <= -5.6e0) goto S60;
00071     if(*ind != 0) goto S30;
00072     if(*x > 100.0e0) goto S70;
00073     if(*x**x > -exparg(&K1)) goto S70;
00074 S30:
00075     t = pow(1.0e0/ *x,2.0);
00076     top = (((r[0]*t+r[1])*t+r[2])*t+r[3])*t+r[4];
00077     bot = (((s[0]*t+s[1])*t+s[2])*t+s[3])*t+1.0e0;
00078     erfc1 = (c-t*top/bot)/ax;
00079 S40:
00080 /*
00081                       FINAL ASSEMBLY
00082 */
00083     if(*ind == 0) goto S50;
00084     if(*x < 0.0e0) erfc1 = 2.0e0*exp(*x**x)-erfc1;
00085     return erfc1;
00086 S50:
00087     w = *x**x;
00088     t = w;
00089     e = w-t;
00090     erfc1 = (0.5e0+(0.5e0-e))*exp(-t)*erfc1;
00091     if(*x < 0.0e0) erfc1 = 2.0e0-erfc1;
00092     return erfc1;
00093 S60:
00094 /*
00095              LIMIT VALUE FOR LARGE NEGATIVE X
00096 */
00097     erfc1 = 2.0e0;
00098     if(*ind != 0) erfc1 = 2.0e0*exp(*x**x);
00099     return erfc1;
00100 S70:
00101 /*
00102              LIMIT VALUE FOR LARGE POSITIVE X
00103                        WHEN IND = 0
00104 */
00105     erfc1 = 0.0e0;
00106     return erfc1;
00107 } /* END */