home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #26 / NN_1992_26.iso / spool / comp / graphics / gnuplot / 314 < prev    next >
Encoding:
Text File  |  1992-11-07  |  2.4 KB  |  72 lines
  1. x-gateway: rodan.UU.NET from info-gnuplot to comp.graphics.gnuplot; Fri, 6 Nov 1992 11:46:24 EST
  2. Message-ID: <199211061349.AA24351@ames.arc.nasa.gov>
  3. Date: Fri, 06 Nov 1992 17:38:00 +0100
  4. From: "MSFD1::MOEHRING" <moehring@msfd1.dnet.gwdg.de>
  5. Newsgroups: comp.graphics.gnuplot
  6. Path: sparky!uunet!wendy-fate.uu.net!info-gnuplot
  7. Sender: info-gnuplot-request@arc.nasa.gov
  8. Subject: (no subject given!)
  9. Lines: 61
  10.  
  11. I have here the version 3.2 of theGNUPLOT program and use it on an MSDOS
  12. computer. I think that the calculation of the complex gamma function is
  13. not too difficult and can easily be included in GNUPLOT. Here is a routine
  14. to compute the gamma-function for complex arguments which is suitable for
  15. inclusion in GNUPLOT.
  16.  
  17.  
  18. /* This approximation is based on Stirlings formula, see e.g.
  19.  M. Abramowitz, I.A. Stegun, Handbook of mathematical functions */
  20.  
  21. f_gamma()
  22.  {
  23.  struct value a;
  24.  static double sti[]={1./12,-1./360,1./1260,-1./1680,5/(66.*90),
  25.               -691/(132.*2730),7/(6.*182)};
  26.  double xr,xi,gr,gi,z1,z2,z3,z4,z5,z6,z7;
  27.  int    n,i;
  28.  char spiegel;
  29.   (void)pop(&a);xr=real(&a);xi=imag(&a);
  30.   spiegel=1;n=0;
  31.   if(xr<0) xr=1-xr; else spiegel=0;
  32.   z1=196-14*xi;
  33.   if((xr-1)*(xr-1)<z1) n=(int)(sqrt(z1)-xr)+2;
  34.   xr+=n;
  35.   if(xr<1) {n=1;xr++;}
  36.   z3=xr*xr+xi*xi;z4=0.5*log(z3);z5=atan(xi/xr);
  37.   gr=(xr-0.5)*z4-xi*z5-xr+0.5*log(2*Pi);
  38.   gi=(xr-0.5)*z5+xi*z4-xi;
  39.   z4=-xi/z3;z3=xr/z3;
  40.   z1=z3*z3-z4*z4;z2=2*z3*z4;
  41.   z5=sti[6];z6=0;
  42.   for(i=5;i>=0;i--)
  43.    {z7=z5*z1-z6*z2;z6=z5*z2+z6*z1;z5=z7+sti[i];}
  44.   gr += z5*z3 - z6*z4;gi+=z5*z4+z6*z3;
  45.   if (gr>88){undefined = TRUE;z5=gr;} else z5=exp(gr);
  46.   gr=z5*cos(gi);gi=z5*sin(gi);
  47.   if(n>0)
  48.    { z1=xr-1;z2=xi;
  49.      for(i=2;i<=n;i++)
  50.       {z3=z1*(xr-i)-z2*xi;z2=z1*xi+z2*(xr-i);z1=z3;}
  51.      z3=z1*z1+z2*z2;z4=(gr*z1+gi*z2)/z3;
  52.      gi=(gi*z1-gr*z2)/z3;gr=z4;
  53.    }
  54.    if(spiegel)
  55.     {z3=exp(-Pi*xi);
  56.      z1=0.5*sin(Pi*(xr-n))*(z3+1/z3);z2=-0.5*cos(Pi*(xr-n))*(z3-1/z3);
  57.      z3=z1*gr-z2*gi;z2=z1*gi+z2*gr;z1=z3;
  58.      if(z3==0.0){undefined = TRUE; z3=1.;} else z3=z1*z1+z2*z2;
  59.      gr=Pi*z1/z3;gi=Pi*z2/z3;
  60.     }
  61.   push( complex(&a,gr,gi));
  62.  }
  63.  
  64. To include this routine in the GNUPLOT program, I have modified the files
  65. "STANDARD.C" and "PLOT.H". In "STANDARD.C" I have removed all the stuff
  66. about the  gamma-function and replaced it by the above function. Then I have
  67. defined GAMMA in "PLOT.H" by:
  68.  
  69. #define GAMMA gamma
  70.  
  71. Then it works fine after compilation with Turbo C.
  72.