home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume18 / planet / part01 / src / rain.c < prev   
Encoding:
C/C++ Source or Header  |  1991-04-10  |  7.0 KB  |  161 lines

  1. /* This program is Copyright (c) 1991 David Allen.  It may be freely
  2.    distributed as long as you leave my name and copyright notice on it.
  3.    I'd really like your comments and feedback; send e-mail to
  4.    allen@viewlogic.com, or send us-mail to David Allen, 10 O'Moore
  5.    Avenue, Maynard, MA 01754. */
  6.  
  7. /* This file contains the routines to compute rainfall. */
  8.  
  9. #include "const.h"
  10. #include "clim.h"
  11.  
  12. /* The input data arrays are l and lm, from main.c, wd, from wind.c,
  13.    and pr, from pressure.c.  Output arrays are rm and rn; fr and fs are
  14.    used as temporary storage. */
  15.  
  16. extern unsigned char l[MAXX][MAXY], lm[MAXX][MAXY], wd[MAXB][MAXX][MAXY];
  17. extern unsigned char pr[MAXB][MAXX][MAXY];
  18. unsigned char rn[MAXB][MAXX][MAXY];
  19. static unsigned char fr[2][MAXX][MAXY], fs[MAXX][MAXY];
  20.  
  21.  
  22. /* The externs below are parameters defined in main.c; the ints are
  23.    parameters defined here. */
  24.  
  25. extern int BSIZE, XSIZE, YSIZE;
  26. int MAXFETCH = 5, RAINCONST = 32, LANDEL = 10, MOUNTDEL = 32, FETCHDEL = 4;
  27. int NRFDEL = 3, HEQDEL = 32, NRHEQDEL = 24, FLANKDEL = -24, PRINTRAIN = 0;
  28.  
  29.  
  30. rainpar (s) char *s; {
  31.    /* This function is called by mainpar() in main.c; it simply tests input
  32.    parameter names to see if they are defined in this file.  Each of the
  33.    above ints are defined in this file.  If the input string matches here,
  34.    the function returns true. */
  35.  
  36.    if      (CMP ("MAXFETCH"))  getlng  (&MAXFETCH,  M_RAIN);
  37.    else if (CMP ("RAINCONST")) getlng  (&RAINCONST, M_RAIN);
  38.    else if (CMP ("LANDEL"))    getlng  (&LANDEL,    M_RAIN);
  39.    else if (CMP ("MOUNTDEL"))  getlng  (&MOUNTDEL,  M_RAIN);
  40.    else if (CMP ("FETCHDEL"))  getlng  (&FETCHDEL,  M_RAIN);
  41.    else if (CMP ("HEQDEL"))    getlng  (&HEQDEL,    M_RAIN);
  42.    else if (CMP ("NRHEQDEL"))  getlng  (&NRHEQDEL,  M_RAIN);
  43.    else if (CMP ("FLANKDEL"))  getlng  (&FLANKDEL,  M_RAIN);
  44.    else if (CMP ("NRFDEL"))    getlng  (&NRFDEL,    M_RAIN);
  45.    else if (CMP ("PRINTRAIN")) getlng  (&PRINTRAIN, M_RAIN);
  46.    else return (0);
  47.    return (1); }
  48.  
  49.  
  50. raincomp () {
  51.    /* This is the main rain computation function.  It calls the functions
  52.    getfetch () and getrain () to do all the work for each buffer, then
  53.    prints out the results if needed. */
  54.  
  55.    register int buf;
  56.  
  57.    for (buf=0; buf<BSIZE; buf++) {
  58.       status (M_RAIN, buf); checkmouse (); 
  59.       getfetch (buf); checkmouse (); getrain (buf); }
  60.  
  61.    if (PRINTRAIN) for (buf=0; buf<BSIZE; buf++)
  62.       putmat ("RAIN", buf, PRINTMODE_SCALE, rn[buf], lm); }
  63.  
  64.  
  65. raindraw (n) int n; { draw (DRAW_GREY, LINE_CORN, rn[n], lm); }
  66.    /* This function calls draw with the right arguments to display rain */
  67.  
  68.  
  69. fetchinc (x, y, dest) int x, y, dest; {
  70.    /* This is the workhorse function for getfetch(), below.  It is called
  71.    several times per square.  It changes x to account for wraparound, so it
  72.    won't work as a macro.  If y is out of range it does nothing, else it
  73.    "marks" the new square in fr[dest] and increments fs to record the number
  74.    of times the square has been marked. */
  75.  
  76.    if (x == -1) x = XSIZE-1; else if (x == XSIZE) x = 0;
  77.    if ((y == -1) || (y == YSIZE)) return (0);
  78.    fr[dest][x][y] = 1; fs[x][y]++; return (0); }
  79.  
  80.  
  81. getfetch (buf) int buf; {
  82.    /* "Fetch" is the term that describes how many squares a given wind line
  83.    travels over water.  It measures how moist the wind is.  The algorithm to
  84.    measure fetch looks like many simultaneous tree walks, where each water
  85.    square is a root square, and every wind edge is a tree edge.  A counter
  86.    for each square determines how many times that square is reached during
  87.    the tree walks; that is the fetch. */
  88.  
  89.    register int i, j, k; int src, dest;
  90.  
  91.    /* Initialize the counter fs to zero.  Array fr, which records the */
  92.    /* list of active edges in the walks, is set so that all ocean squares */
  93.    /* are active.  Also, the result array rn is cleared. */
  94.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
  95.       fr[0][i][j] = l[i][j] ? 0 : 1; fs[i][j] = 0; rn[buf][i][j] = 0; }
  96.  
  97.    /* Each time through the loop, each square is examined.  If it's */
  98.    /* active, disable the mark in the current time step (thus ensuring */
  99.    /* that when the buffers are flipped, the new destination is empty). */
  100.    /* If the square is a mountain, don't pass the mark, but instead add */
  101.    /* some amount to the square -- implementing rain shadows and rainy */
  102.    /* mountain squares.  Finally, for each of the eight cardinal */
  103.    /* directions, if there is wind blowing in that direction, carry a */
  104.    /* marker to that square using fetchinc(), above. */
  105.  
  106.    for (k=0; k<MAXFETCH; k++) {
  107.       src = k % 2; dest = 1 - src;
  108.       for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) if (fr[src][i][j]) {
  109.          fr[src][i][j] = 0;
  110.          if (l[i][j] == 2) rn[buf][i][j] += MOUNTDEL;
  111.          else switch (wd[buf][i][j]) {
  112.             case N|E: fetchinc (i+1, j-1, dest); break;
  113.             case N|W: fetchinc (i-1, j-1, dest); break;
  114.             case S|E: fetchinc (i+1, j+1, dest); break;
  115.             case S|W: fetchinc (i-1, j+1, dest); break;
  116.             case N:   fetchinc (i,   j-1, dest); break;
  117.             case S:   fetchinc (i,   j+1, dest); break;
  118.             case E:   fetchinc (i+1, j,   dest); break;
  119.             case W:   fetchinc (i-1, j,   dest); break; } } } }
  120.  
  121.  
  122. /* This macro is called several times per square by getrain(), below.  It
  123.    simply tests the square for several conditions: if the square is on the
  124.    heat equator, itcz is set to one; if the wind blows south in this square,
  125.    it is on the flank of a circular wind zone (and thus less rainy); the local
  126.    rain sum, x, is increased according to the fetch sum in the square. */
  127. #define RAINTEST(xx,yy) \
  128.    if (pr[buf][xx][yy] == PR_HEQ) itcz = 1; \
  129.    if (wd[buf][xx][yy] & S) flank = 1; \
  130.    x += (fs[xx][yy] * NRFDEL);
  131.  
  132.  
  133. getrain (buf) int buf; {
  134.    /* Once the fetch array is computed, this function looks at each square to
  135.    determine the amount of rainfall there.  The above macro is called five
  136.    times, once for the square and each of its four neighbors; this determines
  137.    whether the square is near the ITCZ or the flank of an air cycle.  The
  138.    sum of fetches for the neighbors is also determined.   Finally, each of the
  139.    factors is weighted and added to the rainfall value:  the local fetch value,
  140.    a land factor, the nearness of the heat equator, and the nearness of a
  141.    flank.  Note that while rn is zeroed in getfetch(), it may be increased by
  142.    rain falling on mountains, so it is nonzero when this function is called. */
  143.  
  144.    register int i, j, x; int itcz, flank;
  145.  
  146.    for (i=0; i<XSIZE; i++) for (j=0; j<YSIZE; j++) {
  147.       flank = 0; itcz = 0; x = rn[buf][i][j];
  148.       if (i < XSIZE-1) { RAINTEST (i+1, j) } else { RAINTEST (0, j) }
  149.       if (i) { RAINTEST (i-1, j) } else { RAINTEST (XSIZE-1, j) }
  150.       if (j < YSIZE-1) { RAINTEST (i, j+1) }
  151.       if (j) { RAINTEST (i, j-1) }
  152.       RAINTEST (i, j);
  153.  
  154.       x += (RAINCONST + FETCHDEL * fs[i][j]);
  155.       if (l[i][j]) x += LANDEL;
  156.       if (pr[buf][i][j] == PR_HEQ) x += HEQDEL;
  157.       if (itcz) x += NRHEQDEL;
  158.       if (flank) x += FLANKDEL;
  159.       if (x < 0) x = 0; if (x> 255) x = 255;
  160.       rn[buf][i][j] = x; } }
  161.