home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / utilities / utilsm / netpbmsca / ppm / c / ppmforge < prev    next >
Encoding:
Text File  |  1993-10-04  |  27.4 KB  |  952 lines
  1. /*
  2.  
  3.         Fractal forgery generator for the PPM toolkit
  4.  
  5.     Originally  designed  and  implemented    in December of 1989 by
  6.     John Walker as a stand-alone program for the  Sun  and    MS-DOS
  7.     under  Turbo C.  Adapted in September of 1991 for use with Jef
  8.         Poskanzer's raster toolkit.
  9.  
  10.     References cited in the comments are:
  11.  
  12.         Foley, J. D., and Van Dam, A., Fundamentals of Interactive
  13.         Computer  Graphics,  Reading,  Massachusetts:  Addison
  14.         Wesley, 1984.
  15.  
  16.         Peitgen, H.-O., and Saupe, D. eds., The Science Of Fractal
  17.         Images, New York: Springer Verlag, 1988.
  18.  
  19.         Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling,
  20.         W. T., Numerical Recipes In C, New Rochelle: Cambridge
  21.         University Press, 1988.
  22.  
  23.     Author:
  24.         John Walker
  25.         Autodesk SA
  26.         Avenue des Champs-Montants 14b
  27.         CH-2074 MARIN
  28.         Switzerland
  29.         Usenet: kelvin@Autodesk.com
  30.         Fax:    038/33 88 15
  31.         Voice:  038/33 76 33
  32.  
  33.     Permission    to  use, copy, modify, and distribute this software and
  34.     its documentation  for  any  purpose  and  without    fee  is  hereby
  35.     granted,  without any conditions or restrictions.  This software is
  36.     provided "as is" without express or implied warranty.
  37.  
  38.                 PLUGWARE!
  39.  
  40.     If you like this kind of stuff, you may also enjoy "James  Gleick's
  41.     Chaos--The  Software"  for  MS-DOS,  available for $59.95 from your
  42.     local software store or directly from Autodesk, Inc., Attn: Science
  43.     Series,  2320  Marinship Way, Sausalito, CA 94965, USA.  Telephone:
  44.     (800) 688-2344 toll-free or, outside the  U.S. (415)  332-2344  Ext
  45.     4886.   Fax: (415) 289-4718.  "Chaos--The Software" includes a more
  46.     comprehensive   fractal    forgery      generator    which    creates
  47.     three-dimensional  landscapes  as  well as clouds and planets, plus
  48.     five more modules which explore other aspects of Chaos.   The  user
  49.     guide  of  more  than  200    pages includes an introduction by James
  50.     Gleick and detailed explanations by Rudy Rucker of the  mathematics
  51.     and algorithms used by each program.
  52.  
  53. */
  54.  
  55. #include <math.h>
  56. #include "ppm.h"
  57.  
  58. #ifndef M_PI
  59. #define M_PI    3.14159265358979323846
  60. #endif
  61. #ifndef M_E
  62. #define M_E    2.7182818284590452354
  63. #endif
  64.  
  65. /* Definitions used to address real and imaginary parts in a two-dimensional
  66.    array of complex numbers as stored by fourn(). */
  67.  
  68. #define Real(v, x, y)  v[1 + (((x) * meshsize) + (y)) * 2]
  69. #define Imag(v, x, y)  v[2 + (((x) * meshsize) + (y)) * 2]
  70.  
  71. /* Co-ordinate indices within arrays. */
  72.  
  73. #define X     0
  74. #define Y     1
  75. #define Z     2
  76.  
  77. /* Definition for obtaining random numbers. */
  78.  
  79. #define nrand 4               /* Gauss() sample count */
  80. #define Cast(low, high) ((low)+(((high)-(low)) * ((random() & 0x7FFF) / arand)))
  81.  
  82. /* Utility definition to get an  array's  element  count  (at  compile
  83.    time).   For  example:  
  84.  
  85.        int  arr[] = {1,2,3,4,5};
  86.        ... 
  87.        printf("%d", ELEMENTS(arr));
  88.  
  89.    would print a five.  ELEMENTS("abc") can also be used to  tell  how
  90.    many  bytes are in a string constant INCLUDING THE TRAILING NULL. */
  91.    
  92. #define ELEMENTS(array) (sizeof(array)/sizeof((array)[0]))
  93.  
  94. /*  Data types    */
  95.  
  96. typedef int Boolean;
  97. #define FALSE 0
  98. #define TRUE 1
  99.  
  100. #define V     (void)
  101.  
  102. /*  Display parameters    */
  103.  
  104. #define SCRSAT     255              /* Screen maximum intensity */
  105.  
  106. /* prototypes */
  107. static void fourn ARGS((float data[], int nn[], int ndim, int isign));
  108. static void initgauss ARGS((unsigned int seed));
  109. static double gauss ARGS((void));
  110. static void spectralsynth ARGS((float **x, unsigned int n, double h));
  111. static void initseed ARGS((void));
  112. static void temprgb ARGS((double temp, double *r, double *g, double *b));
  113. static void etoile ARGS((pixel *pix));
  114. static void genplanet ARGS((float *a, unsigned int n));
  115. static Boolean planet ARGS((void));
  116. /*  Local variables  */
  117.  
  118. static double arand, gaussadd, gaussfac; /* Gaussian random parameters */
  119. static double fracdim;              /* Fractal dimension */
  120. static double powscale;           /* Power law scaling exponent */
  121. static int meshsize = 256;          /* FFT mesh size */
  122. static unsigned int rseed;          /* Current random seed */
  123. static Boolean seedspec = FALSE;      /* Did the user specify a seed ? */
  124. static Boolean clouds = FALSE;          /* Just generate clouds */
  125. static Boolean stars = FALSE;          /* Just generate stars */
  126. static int screenxsize = 256;          /* Screen X size */
  127. static int screenysize = 256;          /* Screen Y size */
  128. static double inclangle, hourangle;   /* Star position relative to planet */
  129. static Boolean inclspec = FALSE;      /* No inclination specified yet */
  130. static Boolean hourspec = FALSE;      /* No hour specified yet */
  131. static double icelevel;           /* Ice cap theshold */
  132. static double glaciers;           /* Glacier level */
  133. static int starfraction;          /* Star fraction */
  134. static int starcolour;              /* Star colour saturation */
  135.  
  136. /*    FOURN  --  Multi-dimensional fast Fourier transform
  137.  
  138.     Called with arguments:
  139.  
  140.        data       A  one-dimensional  array  of  floats  (NOTE!!!    NOT
  141.               DOUBLES!!), indexed from one (NOTE!!!   NOT  ZERO!!),
  142.               containing  pairs of numbers representing the complex
  143.               valued samples.  The Fourier transformed results    are
  144.               returned in the same array.
  145.  
  146.        nn          An  array specifying the edge size in each dimension.
  147.               THIS ARRAY IS INDEXED FROM  ONE,    AND  ALL  THE  EDGE
  148.               SIZES MUST BE POWERS OF TWO!!!
  149.  
  150.        ndim       Number of dimensions of FFT to perform.  Set to 2 for
  151.               two dimensional FFT.
  152.  
  153.        isign      If 1, a Fourier transform is done; if -1 the  inverse
  154.               transformation is performed.
  155.  
  156.         This  function  is essentially as given in Press et al., "Numerical
  157.         Recipes In C", Section 12.11, pp.  467-470.
  158. */
  159.  
  160. static void fourn(data, nn, ndim, isign)
  161.   float data[];
  162.   int nn[], ndim, isign;
  163. {
  164.     register int i1, i2, i3;
  165.     int i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2;
  166.     int ibit, idim, k1, k2, n, nprev, nrem, ntot;
  167.     float tempi, tempr;
  168.     double theta, wi, wpi, wpr, wr, wtemp;
  169.  
  170. #define SWAP(a,b) tempr=(a); (a) = (b); (b) = tempr
  171.  
  172.     ntot = 1;
  173.     for (idim = 1; idim <= ndim; idim++)
  174.     ntot *= nn[idim];
  175.     nprev = 1;
  176.     for (idim = ndim; idim >= 1; idim--) {
  177.     n = nn[idim];
  178.     nrem = ntot / (n * nprev);
  179.     ip1 = nprev << 1;
  180.     ip2 = ip1 * n;
  181.     ip3 = ip2 * nrem;
  182.     i2rev = 1;
  183.     for (i2 = 1; i2 <= ip2; i2 += ip1) {
  184.         if (i2 < i2rev) {
  185.         for (i1 = i2; i1 <= i2 + ip1 - 2; i1 += 2) {
  186.             for (i3 = i1; i3 <= ip3; i3 += ip2) {
  187.             i3rev = i2rev + i3 - i2;
  188.             SWAP(data[i3], data[i3rev]);
  189.             SWAP(data[i3 + 1], data[i3rev + 1]);
  190.             }
  191.         }
  192.         }
  193.         ibit = ip2 >> 1;
  194.         while (ibit >= ip1 && i2rev > ibit) {
  195.         i2rev -= ibit;
  196.         ibit >>= 1;
  197.         }
  198.         i2rev += ibit;
  199.     }
  200.     ifp1 = ip1;
  201.     while (ifp1 < ip2) {
  202.         ifp2 = ifp1 << 1;
  203.         theta = isign * (M_PI * 2) / (ifp2 / ip1);
  204.         wtemp = sin(0.5 * theta);
  205.         wpr = -2.0 * wtemp * wtemp;
  206.         wpi = sin(theta);
  207.         wr = 1.0;
  208.         wi = 0.0;
  209.         for (i3 = 1; i3 <= ifp1; i3 += ip1) {
  210.         for (i1 = i3; i1 <= i3 + ip1 - 2; i1 += 2) {
  211.             for (i2 = i1; i2 <= ip3; i2 += ifp2) {
  212.             k1 = i2;
  213.             k2 = k1 + ifp1;
  214.             tempr = wr * data[k2] - wi * data[k2 + 1];
  215.             tempi = wr * data[k2 + 1] + wi * data[k2];
  216.             data[k2] = data[k1] - tempr;
  217.             data[k2 + 1] = data[k1 + 1] - tempi;
  218.             data[k1] += tempr;
  219.             data[k1 + 1] += tempi;
  220.             }
  221.         }
  222.         wr = (wtemp = wr) * wpr - wi * wpi + wr;
  223.         wi = wi * wpr + wtemp * wpi + wi;
  224.         }
  225.         ifp1 = ifp2;
  226.     }
  227.     nprev *= n;
  228.     }
  229. }
  230. #undef SWAP
  231.  
  232. /*  INITGAUSS  --  Initialise random number generators.  As given in
  233.            Peitgen & Saupe, page 77. */
  234.  
  235. static void initgauss(seed)
  236.   unsigned int seed;
  237. {
  238.     /* Range of random generator */
  239.     arand = pow(2.0, 15.0) - 1.0;
  240.     gaussadd = sqrt(3.0 * nrand);
  241.     gaussfac = 2 * gaussadd / (nrand * arand);
  242.     srandom(seed);
  243. }
  244.  
  245. /*  GAUSS  --  Return a Gaussian random number.  As given in Peitgen
  246.            & Saupe, page 77. */
  247.  
  248. static double gauss()
  249. {
  250.     int i;
  251.     double sum = 0.0;
  252.  
  253.     for (i = 1; i <= nrand; i++) {
  254.     sum += (random() & 0x7FFF);
  255.     }
  256.     return gaussfac * sum - gaussadd;
  257. }
  258.  
  259. /*  SPECTRALSYNTH  --  Spectrally  synthesised    fractal  motion in two
  260.                dimensions.  This algorithm is given under  the
  261.                name   SpectralSynthesisFM2D  on  page  108  of
  262.                Peitgen & Saupe. */
  263.  
  264. static void spectralsynth(x, n, h)
  265.   float **x;
  266.   unsigned int n;
  267.   double h;
  268. {
  269.     unsigned bl;
  270.     int i, j, i0, j0, nsize[3];
  271.     double rad, phase, rcos, rsin;
  272.     float *a;
  273.  
  274.     bl = ((((unsigned long) n) * n) + 1) * 2 * sizeof(float);
  275.     a = (float *) calloc(bl, 1);
  276.     if (a == (float *) 0) {
  277.         pm_error("Cannot allocate %d x %d result array (%ld bytes).",
  278.        n, n, bl);
  279.     }
  280.     *x = a;
  281.  
  282.     for (i = 0; i <= n / 2; i++) {
  283.     for (j = 0; j <= n / 2; j++) {
  284.         phase = 2 * M_PI * ((random() & 0x7FFF) / arand);
  285.         if (i != 0 || j != 0) {
  286.         rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
  287.         } else {
  288.         rad = 0;
  289.         }
  290.         rcos = rad * cos(phase);
  291.         rsin = rad * sin(phase);
  292.         Real(a, i, j) = rcos;
  293.         Imag(a, i, j) = rsin;
  294.         i0 = (i == 0) ? 0 : n - i;
  295.         j0 = (j == 0) ? 0 : n - j;
  296.         Real(a, i0, j0) = rcos;
  297.         Imag(a, i0, j0) = - rsin;
  298.     }
  299.     }
  300.     Imag(a, n / 2, 0) = 0;
  301.     Imag(a, 0, n / 2) = 0;
  302.     Imag(a, n / 2, n / 2) = 0;
  303.     for (i = 1; i <= n / 2 - 1; i++) {
  304.     for (j = 1; j <= n / 2 - 1; j++) {
  305.         phase = 2 * M_PI * ((random() & 0x7FFF) / arand);
  306.         rad = pow((double) (i * i + j * j), -(h + 1) / 2) * gauss();
  307.         rcos = rad * cos(phase);
  308.         rsin = rad * sin(phase);
  309.         Real(a, i, n - j) = rcos;
  310.         Imag(a, i, n - j) = rsin;
  311.         Real(a, n - i, j) = rcos;
  312.         Imag(a, n - i, j) = - rsin;
  313.     }
  314.     }
  315.  
  316.     nsize[0] = 0;
  317.     nsize[1] = nsize[2] = n;          /* Dimension of frequency domain array */
  318.     fourn(a, nsize, 2, -1);          /* Take inverse 2D Fourier transform */
  319. }
  320.  
  321. /*  INITSEED  --  Generate initial random seed, if needed.  */
  322.  
  323. static void initseed()
  324. {
  325.     int i = time((long *) 0) ^ 0xF37C;
  326.     V srandom(i);
  327.     for (i = 0; i < 7; i++)
  328.     V random();
  329.     rseed = random();
  330. }
  331.  
  332. /*  TEMPRGB  --  Calculate the relative R, G, and B components    for  a
  333.          black    body  emitting    light  at a given temperature.
  334.          The Planck radiation equation is solved directly  for
  335.          the R, G, and B wavelengths defined for the CIE  1931
  336.          Standard    Colorimetric    Observer.      The    colour
  337.          temperature is specified in degrees Kelvin. */
  338.  
  339. static void temprgb(temp, r, g, b)
  340.   double temp;
  341.   double *r, *g, *b;
  342. {
  343.     double c1 = 3.7403e10,
  344.        c2 = 14384.0,
  345.        er, eg, eb, es;
  346.  
  347. /* Lambda is the wavelength in microns: 5500 angstroms is 0.55 microns. */
  348.  
  349. #define Planck(lambda)  ((c1 * pow((double) lambda, -5.0)) /  \
  350.              (pow(M_E, c2 / (lambda * temp)) - 1))
  351.  
  352.     er = Planck(0.7000);
  353.     eg = Planck(0.5461);
  354.     eb = Planck(0.4358);
  355. #undef Planck
  356.  
  357.     es = 1.0 / max(er, max(eg, eb));
  358.  
  359.     *r = er * es;
  360.     *g = eg * es;
  361.     *b = eb * es;
  362. }
  363.  
  364. /*  ETOILE  --    Set a pixel in the starry sky.    */
  365.  
  366. static void etoile(pix)
  367.   pixel *pix;
  368. {
  369.     if ((random() % 1000) < starfraction) {
  370. #define StarQuality    0.5          /* Brightness distribution exponent */
  371. #define StarIntensity    8          /* Brightness scale factor */
  372. #define StarTintExp    0.5          /* Tint distribution exponent */
  373.     double v = StarIntensity * pow(1 / (1 - Cast(0, 0.9999)),
  374.                        (double) StarQuality),
  375.            temp,
  376.            r, g, b;
  377.  
  378.     if (v > 255) {
  379.         v = 255;
  380.     }
  381.  
  382.     /* We make a special case for star colour  of zero in order to
  383.        prevent  floating  point  roundoff  which  would  otherwise
  384.        result  in  more  than  256 star colours.  We can guarantee
  385.        that if you specify no star colour, you never get more than
  386.        256 shades in the image. */
  387.  
  388.     if (starcolour == 0) {
  389.         int vi = v;
  390.  
  391.         PPM_ASSIGN(*pix, vi, vi, vi);
  392.     } else {
  393.         temp = 5500 + starcolour *
  394.               pow(1 / (1 - Cast(0, 0.9999)), StarTintExp) *
  395.                   ((random() & 7) ? -1 : 1);
  396.         /* Constrain temperature to a reasonable value: >= 2600K 
  397.            (S Cephei/R Andromedae), <= 28,000 (Spica). */
  398.         temp = max(2600, min(28000, temp));
  399.         temprgb(temp, &r, &g, &b);
  400.         PPM_ASSIGN(*pix, (int) (r * v + 0.499),
  401.                  (int) (g * v + 0.499),
  402.                  (int) (b * v + 0.499));
  403.     }
  404.     } else {
  405.     PPM_ASSIGN(*pix, 0, 0, 0);
  406.     }
  407. }
  408.  
  409. /*  GENPLANET  --  Generate planet from elevation array.  */
  410.  
  411. static void genplanet(a, n)
  412.   float *a;
  413.   unsigned int n;
  414. {
  415.     int i, j;
  416.     unsigned char *cp, *ap;
  417.     double *u, *u1;
  418.     unsigned int *bxf, *bxc;
  419.  
  420. #define RGBQuant    255
  421.     pixel *pixels;              /* Pixel vector */
  422.     pixel *rpix;              /* Current pixel pointer */
  423.  
  424. #define Atthick 1.03              /* Atmosphere thickness as a percentage
  425.                                          of planet's diameter */
  426.     double athfac = sqrt(Atthick * Atthick - 1.0);
  427.     double sunvec[3];
  428.  
  429.     Boolean flipped = FALSE;
  430.     double shang, siang;
  431.  
  432.     ppm_writeppminit(stdout, screenxsize, screenysize,
  433.              (pixval) RGBQuant, FALSE);
  434.  
  435.     if (!stars) {
  436.     u = (double *) malloc((unsigned int) (screenxsize * sizeof(double)));
  437.     u1 = (double *) malloc((unsigned int) (screenxsize * sizeof(double)));
  438.     bxf = (unsigned int *) malloc((unsigned int) screenxsize *
  439.           sizeof(unsigned int));
  440.     bxc = (unsigned int *) malloc((unsigned int) screenxsize *
  441.           sizeof(unsigned int));
  442.  
  443.     if (u == (double *) 0 || u1 == (double *) 0 ||
  444.         bxf == (unsigned int *) 0 || bxc == (unsigned int *) 0) {
  445.             pm_error("Cannot allocate %d element interpolation tables.",
  446.              screenxsize);
  447.     }
  448.  
  449.     /* Compute incident light direction vector. */
  450.  
  451.     shang = hourspec ? hourangle : Cast(0, 2 * M_PI);
  452.     siang = inclspec ? inclangle : Cast(-M_PI * 0.12, M_PI * 0.12);
  453.  
  454.     sunvec[X] = sin(shang) * cos(siang);
  455.     sunvec[Y] = sin(siang);
  456.     sunvec[Z] = cos(shang) * cos(siang);
  457.  
  458.     /* Allow only 25% of random pictures to be crescents */
  459.  
  460.     if (!hourspec && ((random() % 100) < 75)) {
  461.         flipped = sunvec[Z] < 0 ? TRUE : FALSE;
  462.         sunvec[Z] = abs(sunvec[Z]);
  463.     }
  464.         pm_message("%s: -seed %d -dimension %.2f -power %.2f -mesh %d",
  465.             clouds ? "clouds" : "planet",
  466.         rseed, fracdim, powscale, meshsize);
  467.     if (!clouds) {
  468.         pm_message(
  469.                "        -inclination %.0f -hour %d -ice %.2f -glaciers %.2f",
  470.            (siang * (180.0 / M_PI)),
  471.            (int) (((shang * (12.0 / M_PI)) + 12 +
  472.           (flipped ? 12 : 0)) + 0.5) % 24, icelevel, glaciers);
  473.             pm_message("        -stars %d -saturation %d.",
  474.         starfraction, starcolour);
  475.     }
  476.  
  477.     /* Prescale the grid points into intensities. */
  478.  
  479.     cp = (unsigned char *) malloc(n * n);
  480.     if (cp == (unsigned char *) 0)
  481.         return;
  482.     ap = cp;
  483.     for (i = 0; i < n; i++) {
  484.         for (j = 0; j < n; j++) {
  485.         *ap++ = (255.0 * (Real(a, i, j) + 1.0)) / 2.0;
  486.         }
  487.     }
  488.  
  489.     /* Fill the screen from the computed  intensity  grid  by  mapping
  490.        screen  points onto the grid, then calculating each pixel value
  491.        by bilinear interpolation from  the    surrounding  grid  points.
  492.        (N.b. the pictures would undoubtedly look better when generated
  493.        with small grids if    a  more  well-behaved  interpolation  were
  494.        used.)
  495.  
  496.        Before  we get started, precompute the line-level interpolation
  497.            parameters and store them in an array so we don't  have  to  do
  498.        this every time around the inner loop. */
  499.  
  500. #define UPRJ(a,size) ((a)/((size)-1.0))
  501.  
  502.     for (j = 0; j < screenxsize; j++) {
  503.         double bx = (n - 1) * UPRJ(j, screenxsize);
  504.  
  505.         bxf[j] = floor(bx);
  506.         bxc[j] = bxf[j] + 1;
  507.         u[j] = bx - bxf[j];
  508.         u1[j] = 1 - u[j];
  509.     }
  510.     } else {
  511.         pm_message("night: -seed %d -stars %d -saturation %d.",
  512.         rseed, starfraction, starcolour);
  513.     }
  514.  
  515.     pixels = ppm_allocrow(screenxsize);
  516.     for (i = 0; i < screenysize; i++) {
  517.     double t, t1, by, dy, dysq, sqomdysq, icet, svx, svy, svz,
  518.            azimuth;
  519.     int byf, byc, lcos;
  520.  
  521.     if (!stars) {              /* Skip all this setup if just stars */
  522.         by = (n - 1) * UPRJ(i, screenysize);
  523.         dy = 2 * (((screenysize / 2) - i) / ((double) screenysize));
  524.         dysq = dy * dy;
  525.         sqomdysq = sqrt(1.0 - dysq);
  526.         svx = sunvec[X];
  527.         svy = sunvec[Y] * dy;
  528.         svz = sunvec[Z] * sqomdysq;
  529.         byf = floor(by) * n;
  530.         byc = byf + n;
  531.         t = by - floor(by);
  532.         t1 = 1 - t;
  533.     }
  534.  
  535.     if (clouds) {
  536.  
  537.         /* Render the FFT output as clouds. */
  538.  
  539.         for (j = 0; j < screenxsize; j++) {
  540.         double r = t1 * u1[j] * cp[byf + bxf[j]] +
  541.                t  * u1[j] * cp[byc + bxf[j]] +
  542.                t  * u[j]  * cp[byc + bxc[j]] +
  543.                t1 * u[j]  * cp[byf + bxc[j]];
  544.         pixval w = (r > 127.0) ? (RGBQuant * ((r - 127.0) / 128.0)) :
  545.                0;
  546.  
  547.         PPM_ASSIGN(*(pixels + j), w, w, RGBQuant);
  548.         }
  549.     } else if (stars) {
  550.  
  551.         /* Generate a starry sky.  Note  that no FFT is performed;
  552.            the output is  generated  directly  from  a  power  law
  553.            mapping    of  a  pseudorandom sequence into intensities. */
  554.  
  555.         for (j = 0; j < screenxsize; j++) {
  556.         etoile(pixels + j);
  557.         }
  558.     } else {
  559.         for (j = 0; j < screenxsize; j++) {
  560.         PPM_ASSIGN(*(pixels + j), 0, 0, 0);
  561.         }
  562.         azimuth = asin(((((double) i) / (screenysize - 1)) * 2) - 1);
  563.         icet = pow(abs(sin(azimuth)), 1.0 / icelevel) - 0.5;
  564.         lcos = (screenysize / 2) * abs(cos(azimuth));
  565.         rpix = pixels + (screenxsize / 2) - lcos;
  566.         for (j = (screenxsize / 2) - lcos;
  567.          j <= (screenxsize / 2) + lcos; j++) {
  568.         double r = t1 * u1[j] * cp[byf + bxf[j]] +
  569.                t  * u1[j] * cp[byc + bxf[j]] +
  570.                t  * u[j]  * cp[byc + bxc[j]] +
  571.                t1 * u[j]  * cp[byf + bxc[j]],
  572.                ice;
  573.         int ir, ig, ib;
  574.         static unsigned char pgnd[][3] = {
  575.            {206, 205, 0}, {208, 207, 0}, {211, 208, 0},
  576.            {214, 208, 0}, {217, 208, 0}, {220, 208, 0},
  577.            {222, 207, 0}, {225, 205, 0}, {227, 204, 0},
  578.            {229, 202, 0}, {231, 199, 0}, {232, 197, 0},
  579.            {233, 194, 0}, {234, 191, 0}, {234, 188, 0},
  580.            {233, 185, 0}, {232, 183, 0}, {231, 180, 0},
  581.            {229, 178, 0}, {227, 176, 0}, {225, 174, 0},
  582.            {223, 172, 0}, {221, 170, 0}, {219, 168, 0},
  583.            {216, 166, 0}, {214, 164, 0}, {212, 162, 0},
  584.            {210, 161, 0}, {207, 159, 0}, {205, 157, 0},
  585.            {203, 156, 0}, {200, 154, 0}, {198, 152, 0},
  586.            {195, 151, 0}, {193, 149, 0}, {190, 148, 0},
  587.            {188, 147, 0}, {185, 145, 0}, {183, 144, 0},
  588.            {180, 143, 0}, {177, 141, 0}, {175, 140, 0},
  589.            {172, 139, 0}, {169, 138, 0}, {167, 137, 0},
  590.            {164, 136, 0}, {161, 135, 0}, {158, 134, 0},
  591.            {156, 133, 0}, {153, 132, 0}, {150, 132, 0},
  592.            {147, 131, 0}, {145, 130, 0}, {142, 130, 0},
  593.            {139, 129, 0}, {136, 128, 0}, {133, 128, 0},
  594.            {130, 127, 0}, {127, 127, 0}, {125, 127, 0},
  595.            {122, 127, 0}, {119, 127, 0}, {116, 127, 0},
  596.            {113, 127, 0}, {110, 128, 0}, {107, 128, 0},
  597.            {104, 128, 0}, {102, 127, 0}, { 99, 126, 0},
  598.            { 97, 124, 0}, { 95, 122, 0}, { 93, 120, 0},
  599.            { 92, 117, 0}, { 92, 114, 0}, { 92, 111, 0},
  600.            { 93, 108, 0}, { 94, 106, 0}, { 96, 104, 0},
  601.            { 98, 102, 0}, {100, 100, 0}, {103,    99, 0},
  602.            {106,  99, 0}, {109,  99, 0}, {111, 100, 0},
  603.            {114, 101, 0}, {117, 102, 0}, {120, 103, 0},
  604.            {123, 102, 0}, {125, 102, 0}, {128, 100, 0},
  605.            {130,  98, 0}, {132,  96, 0}, {133,    94, 0},
  606.            {134,  91, 0}, {134,  88, 0}, {134,    85, 0},
  607.            {133,  82, 0}, {131,  80, 0}, {129,    78, 0}
  608.         };
  609.  
  610.         if (r >= 128) {
  611.             int ix = ((r - 128) * (ELEMENTS(pgnd) - 1)) / 127;
  612.  
  613.             /* Land area.  Look up colour based on elevation from
  614.                precomputed colour map table. */
  615.  
  616.             ir = pgnd[ix][0];
  617.             ig = pgnd[ix][1];
  618.             ib = pgnd[ix][2];
  619.         } else {
  620.  
  621.             /* Water.  Generate clouds above water based on
  622.                elevation.  */
  623.  
  624.             ir = ig = r > 64 ? (r - 64) * 4 : 0;
  625.             ib = 255;
  626.         }
  627.  
  628.         /* Generate polar ice caps. */
  629.  
  630.         ice = max(0.0, (icet + glaciers * max(-0.5, (r - 128) / 128.0)));
  631.         if  (ice > 0.125) {
  632.             ir = ig = ib = 255;
  633.         }
  634.  
  635.         /* Apply limb darkening by cosine rule. */
  636.  
  637.         {   double dx = 2 * (((screenxsize / 2) - j) /
  638.                 ((double) screenysize)),
  639.                dxsq = dx * dx,
  640.                ds, di, inx;
  641.             double dsq, dsat;
  642.             di = svx * dx + svy + svz * sqrt(1.0 - dxsq);
  643. #define         PlanetAmbient  0.05
  644.             if (di < 0)
  645.             di = 0;
  646.             di = min(1.0, di);
  647.  
  648.             ds = sqrt(dxsq + dysq);
  649.             ds = min(1.0, ds);
  650.  
  651.             /* Calculate  atmospheric absorption  based on the
  652.                thickness of atmosphere traversed by  light  on
  653.                its way to the surface. */
  654.  
  655. #define         AtSatFac 1.0
  656.             dsq = ds * ds;
  657.             dsat = AtSatFac * ((sqrt(Atthick * Atthick - dsq) -
  658.                 sqrt(1.0 * 1.0 - dsq)) / athfac);
  659. #define         AtSat(x,y) x = ((x)*(1.0-dsat))+(y)*dsat
  660.             AtSat(ir, 127);
  661.             AtSat(ig, 127);
  662.             AtSat(ib, 255);
  663.  
  664.             inx = PlanetAmbient + (1.0 - PlanetAmbient) * di;
  665.             ir *= inx;
  666.             ig *= inx;
  667.             ib *= inx;
  668.         }
  669.  
  670.         PPM_ASSIGN(*rpix, ir, ig, ib);
  671.         rpix++;
  672.         }
  673.  
  674.         /* Left stars */
  675.  
  676. #define StarClose    2
  677.         for (j = 0; j < (screenxsize / 2) - (lcos + StarClose); j++) {
  678.         etoile(pixels + j);
  679.         }
  680.  
  681.         /* Right stars */
  682.  
  683.         for (j = (screenxsize / 2) + (lcos + StarClose);
  684.          j < screenxsize; j++) {
  685.         etoile(pixels + j);
  686.         }
  687.     }
  688.     ppm_writeppmrow(stdout, pixels, screenxsize, RGBQuant, FALSE);
  689.     }
  690.     pm_close(stdout);
  691.  
  692.     ppm_freerow(pixels);
  693.     if (!stars) {
  694.     free((char *) cp);
  695.     free((char *) u);
  696.     free((char *) u1);
  697.     free((char *) bxf);
  698.     free((char *) bxc);
  699.     }
  700. }
  701.  
  702. /*  PLANET  --    Make a planet.    */
  703.  
  704. static Boolean planet()
  705. {
  706.     float *a = (float *) 0;
  707.     int i, j;
  708. #ifdef VMS
  709.     double rmin = HUGE_VAL, rmax = -HUGE_VAL, rmean, rrange;
  710. #else
  711.     double rmin = 1e50, rmax = -1e50, rmean, rrange;
  712. #endif
  713.  
  714.     if (!seedspec) {
  715.     initseed();
  716.     }
  717.     initgauss(rseed);
  718.  
  719.     if (!stars) {
  720.  
  721.     spectralsynth(&a, meshsize, 3.0 - fracdim);
  722.     if (a == (float *) 0) {
  723.         return FALSE;
  724.     }
  725.  
  726.     /* Apply power law scaling if non-unity scale is requested. */
  727.  
  728.     if (powscale != 1.0) {
  729.         for (i = 0; i < meshsize; i++) {
  730.         for (j = 0; j < meshsize; j++) {
  731.            double r = Real(a, i, j);
  732.  
  733.             if (r > 0) {
  734.             Real(a, i, j) = pow(r, powscale);
  735.             }
  736.         }
  737.         }
  738.     }
  739.  
  740.     /* Compute extrema for autoscaling. */
  741.  
  742.     for (i = 0; i < meshsize; i++) {
  743.         for (j = 0; j < meshsize; j++) {
  744.         double r = Real(a, i, j);
  745.  
  746.         rmin = min(rmin, r);
  747.         rmax = max(rmax, r);
  748.         }
  749.     }
  750.     rmean = (rmin + rmax) / 2;
  751.     rrange = (rmax - rmin) / 2;
  752.     for (i = 0; i < meshsize; i++) {
  753.         for (j = 0; j < meshsize; j++) {
  754.         Real(a, i, j) = (Real(a, i, j) - rmean) / rrange;
  755.         }
  756.     }
  757.     }
  758.     genplanet(a, meshsize);
  759.     if (a != (float *) 0) {
  760.     free((char *) a);
  761.     }
  762.     return TRUE;
  763. }
  764.  
  765. /*  MAIN  --  Main program.  */
  766.  
  767. int main(argc, argv)
  768.   int argc;
  769.   char *argv[];
  770. {
  771.     int i;
  772.     char *usage = "\n\
  773.       [-width|-xsize <x>] [-height|-ysize <y>] [-mesh <n>]\n\
  774.       [-clouds] [-dimension <f>] [-power <f>] [-seed <n>]\n\
  775.       [-hour <f>] [-inclination|-tilt <f>] [-ice <f>] [-glaciers <f>]\n\
  776.       [-night] [-stars <n>] [-saturation <n>]";
  777.     Boolean dimspec = FALSE, meshspec = FALSE, powerspec = FALSE,
  778.         widspec = FALSE, hgtspec = FALSE, icespec = FALSE,
  779.         glacspec = FALSE, starspec = FALSE, starcspec = FALSE;
  780.  
  781.  
  782.     ppm_init(&argc, argv);
  783.     i = 1;
  784.     while ((i < argc) && (argv[i][0] == '-') && (argv[i][1] != '\0')) {
  785.  
  786.         if (pm_keymatch(argv[i], "-clouds", 2)) {
  787.         clouds = TRUE;
  788.         } else if (pm_keymatch(argv[i], "-night", 2)) {
  789.         stars = TRUE;
  790.         } else if (pm_keymatch(argv[i], "-dimension", 2)) {
  791.         if (dimspec) {
  792.                 pm_error("already specified a dimension");
  793.         }
  794.         i++;
  795.             if ((i == argc) || (sscanf(argv[i], "%lf", &fracdim)  != 1))
  796.         pm_usage(usage);
  797.         if (fracdim <= 0.0) {
  798.                 pm_error("fractal dimension must be greater than 0");
  799.         }
  800.         dimspec = TRUE;
  801.         } else if (pm_keymatch(argv[i], "-hour", 3)) {
  802.         if (hourspec) {
  803.                 pm_error("already specified an hour");
  804.         }
  805.         i++;
  806.             if ((i == argc) || (sscanf(argv[i], "%lf", &hourangle) != 1))
  807.         pm_usage(usage);
  808.          hourangle = (M_PI / 12.0) * (hourangle + 12.0);
  809.          hourspec = TRUE;
  810.         } else if (pm_keymatch(argv[i], "-inclination", 3) ||
  811.                    pm_keymatch(argv[i], "-tilt", 2)) {
  812.         if (inclspec) {
  813.                 pm_error("already specified an inclination/tilt");
  814.         }
  815.         i++;
  816.             if ((i == argc) || (sscanf(argv[i], "%lf", &inclangle) != 1))
  817.         pm_usage(usage);
  818.         inclangle = (M_PI / 180.0) * inclangle;
  819.         inclspec = TRUE;
  820.         } else if (pm_keymatch(argv[i], "-mesh", 2)) {
  821.         unsigned int j;
  822.  
  823.         if (meshspec) {
  824.                 pm_error("already specified a mesh size");
  825.         }
  826.         i++;
  827.             if ((i == argc) || (sscanf(argv[i], "%d", &meshsize) != 1))
  828.         pm_usage(usage);
  829.  
  830.         /* Force FFT mesh to the next larger power of 2. */
  831.  
  832.         for (j = meshsize; (j & 1) == 0; j >>= 1) ;
  833.  
  834.         if (j != 1) {
  835.         for (j = 2; j < meshsize; j <<= 1) ;
  836.         meshsize = j;
  837.         }
  838.         meshspec = TRUE;
  839.         } else if (pm_keymatch(argv[i], "-power", 2)) {
  840.         if (powerspec) {
  841.                 pm_error("already specified a power factor");
  842.         }
  843.         i++;
  844.             if ((i == argc) || (sscanf(argv[i], "%lf", &powscale) != 1))
  845.         pm_usage(usage);
  846.         if (powscale <= 0.0) {
  847.                 pm_error("power factor must be greater than 0");
  848.         }
  849.         powerspec = TRUE;
  850.         } else if (pm_keymatch(argv[i], "-ice", 3)) {
  851.         if (icespec) {
  852.                 pm_error("already specified ice cap level");
  853.         }
  854.         i++;
  855.             if ((i == argc) || (sscanf(argv[i], "%lf", &icelevel) != 1))
  856.         pm_usage(usage);
  857.         if (icelevel <= 0.0) {
  858.                 pm_error("ice cap level must be greater than 0");
  859.         }
  860.         icespec = TRUE;
  861.         } else if (pm_keymatch(argv[i], "-glaciers", 2)) {
  862.         if (glacspec) {
  863.                 pm_error("already specified glacier level");
  864.         }
  865.         i++;
  866.             if ((i == argc) || (sscanf(argv[i], "%lf", &glaciers) != 1))
  867.         pm_usage(usage);
  868.         if (glaciers <= 0.0) {
  869.                 pm_error("glacier level must be greater than 0");
  870.         }
  871.         glacspec = TRUE;
  872.         } else if (pm_keymatch(argv[i], "-stars", 3)) {
  873.         if (starspec) {
  874.                 pm_error("already specified a star fraction");
  875.         }
  876.         i++;
  877.             if ((i == argc) || (sscanf(argv[i], "%d", &starfraction) != 1))
  878.         pm_usage(usage);
  879.         starspec = TRUE;
  880.         } else if (pm_keymatch(argv[i], "-saturation", 3)) {
  881.         if (starcspec) {
  882.                 pm_error("already specified a star colour saturation");
  883.         }
  884.         i++;
  885.             if ((i == argc) || (sscanf(argv[i], "%d", &starcolour) != 1))
  886.         pm_usage(usage);
  887.         starcspec = TRUE;
  888.         } else if (pm_keymatch(argv[i], "-seed", 3)) {
  889.         if (seedspec) {
  890.                 pm_error("already specified a random seed");
  891.         }
  892.         i++;
  893.             if ((i == argc) || (sscanf(argv[i], "%d", &rseed) != 1))
  894.         pm_usage(usage);
  895.         seedspec = TRUE;
  896.         } else if (pm_keymatch(argv[i], "-xsize", 2) ||
  897.                    pm_keymatch(argv[i], "-width", 2)) {
  898.         if (widspec) {
  899.                 pm_error("already specified a width/xsize");
  900.         }
  901.         i++;
  902.             if ((i == argc) || (sscanf(argv[i], "%d", &screenxsize) != 1))
  903.         pm_usage(usage);
  904.         widspec = TRUE;
  905.         } else if (pm_keymatch(argv[i], "-ysize", 2) ||
  906.                    pm_keymatch(argv[i], "-height", 3)) {
  907.         if (hgtspec) {
  908.                 pm_error("already specified a height/ysize");
  909.         }
  910.         i++;
  911.             if ((i == argc) || (sscanf(argv[i], "%d", &screenysize) != 1))
  912.         pm_usage(usage);
  913.         hgtspec = TRUE;
  914.     } else {
  915.         pm_usage(usage);
  916.     }
  917.     i++;
  918.     }
  919.  
  920.     /* Set defaults when explicit specifications were not given.
  921.  
  922.        The  default  fractal  dimension  and  power  scale depend upon
  923.        whether we're generating a planet or clouds. */
  924.  
  925.     if (!dimspec) {
  926.     fracdim = clouds ? 2.15 : 2.4;
  927.     }
  928.     if (!powerspec) {
  929.     powscale = clouds ? 0.75 : 1.2;
  930.     }
  931.     if (!icespec) {
  932.     icelevel = 0.4;
  933.     }
  934.     if (!glacspec) {
  935.     glaciers = 0.75;
  936.     }
  937.     if (!starspec) {
  938.     starfraction = 100;
  939.     }
  940.     if (!starcspec) {
  941.     starcolour = 125;
  942.     }
  943.  
  944.     /* Force  screen to be at least  as wide as it is high.  Long,
  945.        skinny screens  cause  crashes  because    picture  width    is
  946.        calculated based on height.  */
  947.  
  948.     screenxsize = max(screenysize, screenxsize);
  949.     screenxsize = (screenxsize + 1) & (~1);
  950.     exit(planet() ? 0 : 1);
  951. }
  952.