archive/o/standard.c

   1 /************************************************************************
   2  *
   3  *  Code is modified from the code for Gnuplot by Zou Maorong
   4  */
   5 /*
   6  *  G N U P L O T  --  standard.c
   7  *  Copyright (C) 1986, 1987  Thomas Williams, Colin Kelley
   8  *  You may use this code as you wish if credit is given and this message
   9  *  is retained.
  10  */
  11 /****************************************************************************/
  12
  13 #include <math.h>
  14 #include <stdio.h>
  15 #include "plot.h"
  16
  17 extern BOOLEAN undefined;
  18 extern int     errno;
  19 extern struct  value stack[STACK_DEPTH];
  20 extern int     s_p;
  21
  22 struct value *pop(), *complex(), *integer();
  23 double magnitude(), angle(), real(), imag();
  24 long   lrand48();
  25 void   srand48();
  26
  27 /***************************************************************************/
  28
  29 f_real()
  30 {
  31 struct value a;
  32         push( complex(&a,real(pop(&a)), 0.0) );
  33 }
  34
  35 f_imag()
  36 {
  37 struct value a;
  38         push( complex(&a,imag(pop(&a)), 0.0) );
  39 }
  40
  41 f_arg()
  42 {
  43 struct value a;
  44         push( complex(&a,angle(pop(&a)), 0.0) );
  45 }
  46
  47 f_conjg()
  48 {
  49 struct value a;
  50         (void) pop(&a);
  51         push( complex(&a,real(&a),-imag(&a) ));
  52 }
  53
  54 f_sin()
  55 {
  56 struct value a;
  57         (void) pop(&a);
  58         push( complex(&a,sin(real(&a))*cosh(imag(&a)), cos(real(&a))*sinh(imag(&a))) );
  59 }
  60
  61 f_cos()
  62 {
  63 struct value a;
  64         (void) pop(&a);
  65         push( complex(&a,cos(real(&a))*cosh(imag(&a)), -sin(real(&a))*sinh(imag(&a))));
  66 }
  67
  68 f_tan()
  69 {
  70 struct value a;
  71 register double den;
  72         (void) pop(&a);
  73         if (imag(&a) == 0.0)
  74                 push( complex(&a,tan(real(&a)),0.0) );
  75         else {
  76                 den = cos(2*real(&a))+cosh(2*imag(&a));
  77                 if (den == 0.0) {
  78                         undefined = TRUE;
  79                         push( &a );
  80                 }
  81                 else
  82                         push( complex(&a,sin(2*real(&a))/den, sinh(2*imag(&a))/den) );
  83         }
  84 }
  85
  86 f_asin()
  87 {
  88 struct value a;
  89 register double alpha, beta, x, y;
  90         (void) pop(&a);
  91         x = real(&a); y = imag(&a);
  92         if (y == 0.0) {
  93                 if (fabs(x) > 1.0) {
  94                         undefined = TRUE;
  95                         push(complex(&a,0.0, 0.0));
  96                 } else
  97                         push( complex(&a,asin(x),0.0) );
  98         } else {
  99                 beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
 100                 alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
 101                 push( complex(&a,asin(beta), log(alpha + sqrt(alpha*alpha-1))) );
 102         }
 103 }
 104
 105 f_acos()
 106 {
 107 struct value a;
 108 register double alpha, beta, x, y;
 109         (void) pop(&a);
 110         x = real(&a); y = imag(&a);
 111         if (y == 0.0) {
 112                 if (fabs(x) > 1.0) {
 113                         undefined = TRUE;
 114                         push(complex(&a,0.0, 0.0));
 115                 } else
 116                         push( complex(&a,acos(x),0.0) );
 117         } else {
 118                 alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
 119                 beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
 120                 push( complex(&a,acos(beta), log(alpha + sqrt(alpha*alpha-1))) );
 121         }
 122 }
 123
 124 f_atan()
 125 {
 126 struct value a;
 127 register double x, y;
 128         (void) pop(&a);
 129         x = real(&a); y = imag(&a);
 130         if (y == 0.0)
 131                 push( complex(&a,atan(x), 0.0) );
 132         else if (x == 0.0 && fabs(y) == 1.0) {
 133                 undefined = TRUE;
 134                 push(complex(&a,0.0, 0.0));
 135         } else
 136                 push( complex(&a,atan(2*x/(1-x*x-y*y)),
 137                         log((x*x+(y+1)*(y+1))/(x*x+(y-1)*(y-1)))/4) );
 138 }
 139
 140 f_sinh()
 141 {
 142 struct value a;
 143         (void) pop(&a);
 144         push( complex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
 145 }
 146
 147 f_cosh()
 148 {
 149 struct value a;
 150         (void) pop(&a);
 151         push( complex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
 152 }
 153
 154 f_tanh()
 155 {
 156 struct value a;
 157 register double den;
 158         (void) pop(&a);
 159         den = cosh(2*real(&a)) + cos(2*imag(&a));
 160         push( complex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
 161 }
 162
 163 f_int()
 164 {
 165 struct value a;
 166         push( integer(&a,(int)real(pop(&a))) );
 167 }
 168
 169
 170 f_abs()
 171 {
 172 struct value a;
 173         (void) pop(&a);
 174         switch (a.type) {
 175                 case INT:
 176                         push( integer(&a,abs(a.v.int_val)) );
 177                         break;
 178                 case CMPLX:
 179                         push( complex(&a,magnitude(&a), 0.0) );
 180         }
 181 }
 182
 183 f_sgn()
 184 {
 185 struct value a;
 186         (void) pop(&a);
 187         switch(a.type) {
 188                 case INT:
 189                         push( integer(&a,(a.v.int_val > 0) ? 1 :
 190                                         (a.v.int_val < 0) ? -1 : 0) );
 191                         break;
 192                 case CMPLX:
 193                         push( integer(&a,(a.v.cmplx_val.real > 0.0) ? 1 :
 194                                         (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
 195                         break;
 196         }
 197 }
 198
 199
 200 f_sqrt()
 201 {
 202 struct value a;
 203 register double mag, ang;
 204         (void) pop(&a);
 205         mag = sqrt(magnitude(&a));
 206         if (imag(&a) == 0.0 && real(&a) < 0.0)
 207                 push( complex(&a,0.0,mag) );
 208         else
 209         {
 210                 if ( (ang = angle(&a)) < 0.0)
 211                         ang += 2*Pi;
 212                 ang /= 2;
 213                 push( complex(&a,mag*cos(ang), mag*sin(ang)) );
 214         }
 215 }
 216
 217
 218 f_exp()
 219 {
 220 struct value a;
 221 register double mag, ang;
 222         (void) pop(&a);
 223         mag = exp(real(&a));
 224         ang = imag(&a);
 225         push( complex(&a,mag*cos(ang), mag*sin(ang)) );
 226 }
 227
 228
 229 f_log10()
 230 {
 231 struct value a;
 232 register double l10;;
 233         (void) pop(&a);
 234         l10 = log(10.0);        /***** replace with a constant! ******/
 235         push( complex(&a,log(magnitude(&a))/l10, angle(&a)/l10) );
 236 }
 237
 238
 239 f_log()
 240 {
 241 struct value a;
 242         (void) pop(&a);
 243         push( complex(&a,log(magnitude(&a)), angle(&a)) );
 244 }
 245
 246
 247 f_besj0()       /* j0(a) = sin(a)/a */
 248 {
 249 struct value a;
 250         a = top_of_stack;
 251         f_sin();
 252         push(&a);
 253         f_div();
 254 }
 255
 256
 257 f_besj1()       /* j1(a) = sin(a)/(a**2) - cos(a)/a */
 258 {
 259 struct value a;
 260         a = top_of_stack;
 261         f_sin();
 262         push(&a);
 263         push(&a);
 264         f_mult();
 265         f_div();
 266         push(&a);
 267         f_cos();
 268         push(&a);
 269         f_div();
 270         f_minus();
 271 }
 272
 273
 274 f_besy0()       /* y0(a) = -cos(a)/a */
 275 {
 276 struct value a;
 277         a = top_of_stack;
 278         f_cos();
 279         push(&a);
 280         f_div();
 281         f_uminus();
 282 }
 283
 284
 285 f_besy1()       /* y1(a) = -cos(a)/(a**2) - sin(a)/a */
 286 {
 287 struct value a;
 288
 289         a = top_of_stack;
 290         f_cos();
 291         push(&a);
 292         push(&a);
 293         f_mult();
 294         f_div();
 295         push(&a);
 296         f_sin();
 297         push(&a);
 298         f_div();
 299         f_plus();
 300         f_uminus();
 301 }
 302
 303
 304 f_floor()
 305 {
 306 struct value a;
 307
 308         (void) pop(&a);
 309         switch (a.type) {
 310                 case INT:
 311                         push( integer(&a,(int)floor((double)a.v.int_val)));
 312                         break;
 313                 case CMPLX:
 314                         push( complex(&a,floor(a.v.cmplx_val.real),
 315                                 floor(a.v.cmplx_val.imag)) );
 316         }
 317 }
 318
 319 f_rand()
 320 {
 321   struct value a;
 322   push( integer(&a,(int)lrand48()));
 323 }
 324
 325 f_srand()
 326 {
 327   struct value a;
 328   (void) pop(&a);
 329   srand48( (long) a.v.int_val);
 330   push( &a);
 331 }
 332
 333
 334 f_ceil()
 335 {
 336 struct value a;
 337
 338         (void) pop(&a);
 339         switch (a.type) {
 340                 case INT:
 341                         push( integer(&a,(int)ceil((double)a.v.int_val)));
 342                         break;
 343                 case CMPLX:
 344                         push( complex(&a,ceil(a.v.cmplx_val.real), ceil(a.v.cmplx_val.imag)) );
 345         }
 346 }
 347
 348 /*
 349 f_gamma()
 350 {
 351 extern int signgam;
 352 register double y;
 353 struct value a;
 354
 355         y = gamma(real(pop(&a)));
 356         if (y > 88.0) {
 357                 undefined = TRUE;
 358                 push( integer(&a,0) );
 359         }
 360         else
 361           push( complex(&a,signgam * exp(y),0.0) );
 362 }
 363
 364 */
 365 /****************************************************************************/