home *** CD-ROM | disk | FTP | other *** search
/ BCI NET / BCI NET Dec 94.iso / archives / programming / source / thesource6.dms / thesource6.adf / Source / Misc / ReactDiff.lha / ReactionDiffusion / spots.c < prev   
Encoding:
C/C++ Source or Header  |  1993-01-29  |  10.3 KB  |  526 lines
  1. /*
  2.  
  3. Make spots and stripes with reaction-diffusion.
  4.  
  5. The spot-formation system is described in the article:
  6.  
  7.   "A Model for Generating Aspects of Zebra and Other Mammailian
  8.    Coat Patterns"
  9.   Jonathan B. L. Bard
  10.   Journal of Theoretical Biology, Vol. 93, No. 2, pp. 363-385
  11.   (November 1981)
  12.  
  13. The stripe-formation system is described in the book:
  14.  
  15.   Models of Biological Pattern Formation
  16.   Hans Meinhardt
  17.   Academic Press, 1982
  18.  
  19.  
  20. Permission is granted to modify and/or distribute this program so long
  21. as the program is distributed free of charge and this header is retained
  22. as part of the program.
  23.  
  24. Copyright (c) Greg Turk, 1991
  25.  
  26. */
  27.  
  28. #include <stdio.h>
  29. #include <math.h>
  30.  
  31. extern double atof();
  32. extern double drand48();
  33.  
  34. float frand();
  35.  
  36. /* screen stuff */
  37.  
  38. int xsize = 30;
  39. int ysize = 30;
  40. int psize = 4;
  41.  
  42. /* simulation variables */
  43.  
  44. int interval = 200;
  45. int value_switch = 1;
  46. int stripe_flag = 0;
  47.  
  48. #define  STRIPES  1
  49. #define  SPOTS    2
  50.  
  51. #define MAX 200
  52.  
  53. float a[MAX][MAX];
  54. float b[MAX][MAX];
  55. float c[MAX][MAX];
  56. float d[MAX][MAX];
  57. float e[MAX][MAX];
  58.  
  59. float da[MAX][MAX];
  60. float db[MAX][MAX];
  61. float dc[MAX][MAX];
  62. float dd[MAX][MAX];
  63. float de[MAX][MAX];
  64.  
  65. float ai[MAX][MAX];
  66.  
  67. float p1,p2,p3;
  68.  
  69. float diff1,diff2;
  70.  
  71. float arand;
  72. float a_steady;
  73. float b_steady;
  74.  
  75. float beta_init;
  76. float beta_rand;
  77.  
  78. float speed = 1.0;
  79.  
  80. int sim = 1;
  81.  
  82.  
  83. /******************************************************************************
  84. Main routine.
  85. ******************************************************************************/
  86.  
  87. main (argc,argv)
  88.   int  argc;
  89.   char *argv[];
  90. {
  91.   char *s;
  92.  
  93.   /* parse the command line options */
  94.  
  95.   while(--argc >0 && (*++argv)[0]=='-') {
  96.     for(s = argv[0]+1; *s; s++)
  97.       switch(*s) {
  98.     case 'n':
  99.       xsize = atoi (*++argv);
  100.       ysize = atoi (*++argv);
  101.           if (xsize > MAX || ysize > MAX)
  102.             printf ("oops, too large a screen size\n");
  103.       argc -= 2;
  104.       break;
  105.     case 'x':
  106.       stripe_flag = 1;
  107.       break;
  108.     default:
  109.       break;
  110.       }
  111.   }
  112.  
  113.   /* setup graphics */
  114.  
  115.   init_graphics (psize * xsize, psize * ysize, 40);
  116.   set_pixel_size (psize);
  117.  
  118.   /* make spots or stripes */
  119.  
  120.   if (stripe_flag)
  121.     do_stripes();
  122.   else
  123.     do_spots();
  124. }
  125.  
  126.  
  127. /******************************************************************************
  128. Run Meinhardt's stripe-formation system.
  129. ******************************************************************************/
  130.  
  131. do_stripes()
  132. {
  133.   p1 = 0.04;
  134.   p2 = 0.06;
  135.   p3 = 0.04;
  136.  
  137.   diff1 = 0.009;
  138.   diff2 = 0.2;
  139.  
  140.   arand = 0.02;
  141.  
  142.   sim = STRIPES;
  143.   value_switch = 1;
  144.   compute();
  145. }
  146.  
  147.  
  148. /******************************************************************************
  149. Run Turing reaction-diffusion system.
  150. ******************************************************************************/
  151.  
  152. do_spots()
  153. {
  154.   beta_init = 12;
  155.   beta_rand = 0.1;
  156.  
  157.   a_steady = 4;
  158.   b_steady = 4;
  159.  
  160.   diff1 = 0.25;
  161.   diff2 = 0.0625;
  162.  
  163.   p1 = 0.2;
  164.   p2 = 0.0;
  165.   p3 = 0.0;
  166.  
  167.   sim = SPOTS;
  168.   value_switch = 2;
  169.   compute();
  170. }
  171.  
  172.  
  173. /******************************************************************************
  174. Diffuse and react.
  175. ******************************************************************************/
  176.  
  177. compute()
  178. {
  179.   int k;
  180.   int iterations = 99999999;
  181.  
  182.   /* calculate semistable equilibria */
  183.  
  184.   semi_equilibria();
  185.  
  186.   /* start things diffusing */
  187.  
  188.   for (k = 0; k < iterations; k++) {
  189.  
  190.     if (k % interval == 0) {
  191.       printf ("iteration %d\n", k);
  192.       switch (value_switch) {
  193.     case 1:
  194.       show(a);
  195.       break;
  196.     case 2:
  197.       show(b);
  198.       break;
  199.     case 3:
  200.       show(c);
  201.       break;
  202.     case 4:
  203.       show(d);
  204.       break;
  205.     case 5:
  206.       show(e);
  207.       break;
  208.     default:
  209.       printf ("bad switch in compute: %d\n", value_switch);
  210.       break;
  211.       }
  212.     }
  213.  
  214.     /* perform reaction and diffusion */
  215.  
  216.     switch (sim) {
  217.       case STRIPES:
  218.     multiplicative_help();
  219.     break;
  220.       case SPOTS:
  221.     turing();
  222.     break;
  223.       default: break;
  224.     }
  225.   }
  226. }
  227.  
  228.  
  229. /******************************************************************************
  230. Create stripes with what Hans Meinhardt calls a two-species balance.
  231. ******************************************************************************/
  232.  
  233. multiplicative_help()
  234. {
  235.   int i,j;
  236.   int iprev,inext,jprev,jnext;
  237.   float aval,bval,cval,dval,eval;
  238.   float ka,kc,kd;
  239.   float temp1,temp2;
  240.   float dda,ddb;
  241.   float ddd,dde;
  242.  
  243.   /* compute change in each cell */
  244.  
  245.   for (i = 0; i < xsize; i++) {
  246.  
  247.     ka = -p1 - 4 * diff1;
  248.     kc = -p2;
  249.     kd = -p3 - 4 * diff2;
  250.  
  251.     iprev = (i + xsize - 1) % xsize;
  252.     inext = (i + 1) % xsize;
  253.  
  254.     for (j = 0; j < ysize; j++) {
  255.  
  256.       jprev = (j + ysize - 1) % ysize;
  257.       jnext = (j + 1) % ysize;
  258.  
  259.       aval = a[i][j];
  260.       bval = b[i][j];
  261.       cval = c[i][j];
  262.       dval = d[i][j];
  263.       eval = e[i][j];
  264.  
  265.       temp1 = 0.01 * aval * aval * eval * ai[i][j];
  266.       temp2 = 0.01 * bval * bval * dval;
  267.  
  268.       dda = a[i][jprev] + a[i][jnext] + a[iprev][j] + a[inext][j];
  269.       ddb = b[i][jprev] + b[i][jnext] + b[iprev][j] + b[inext][j];
  270.       ddd = d[i][jprev] + d[i][jnext] + d[iprev][j] + d[inext][j];
  271.       dde = e[i][jprev] + e[i][jnext] + e[iprev][j] + e[inext][j];
  272.  
  273.       da[i][j] = aval * ka + diff1 * dda + temp1 / cval;
  274.       db[i][j] = bval * ka + diff1 * ddb + temp2 / cval;
  275.       dc[i][j] = cval * kc + temp1 + temp2;
  276.       dd[i][j] = dval * kd + diff2 * ddd + p3 * aval;
  277.       de[i][j] = eval * kd + diff2 * dde + p3 * bval;
  278.     }
  279.   }
  280.  
  281.   /* affect change */
  282.  
  283.   for (i = 0; i < xsize; i++)
  284.     for (j = 0; j < ysize; j++) {
  285.       a[i][j] += (speed * da[i][j]);
  286.       b[i][j] += (speed * db[i][j]);
  287.       c[i][j] += (speed * dc[i][j]);
  288.       d[i][j] += (speed * dd[i][j]);
  289.       e[i][j] += (speed * de[i][j]);
  290.     }
  291. }
  292.  
  293.  
  294. /******************************************************************************
  295. Turing's reaction-diffusion equations.
  296. ******************************************************************************/
  297.  
  298. turing()
  299. {
  300.   int i,j;
  301.   int iprev,inext,jprev,jnext;
  302.   float aval,bval;
  303.   float ka;
  304.   float dda,ddb;
  305.   float Diff1,Diff2;
  306.  
  307.   Diff1 = diff1 / 2.0;
  308.   Diff2 = diff2 / 2.0;
  309.   ka = p1 / 16.0;
  310.  
  311.   /* compute change in each cell */
  312.  
  313.   for (i = 0; i < xsize; i++) {
  314.  
  315.     iprev = (i + xsize - 1) % xsize;
  316.     inext = (i + 1) % xsize;
  317.  
  318.     for (j = 0; j < ysize; j++) {
  319.  
  320.       jprev = (j + ysize - 1) % ysize;
  321.       jnext = (j + 1) % ysize;
  322.  
  323.       aval = a[i][j];
  324.       bval = b[i][j];
  325.  
  326.       dda = a[i][jprev] + a[i][jnext] + a[iprev][j] + a[inext][j] - 4 * aval;
  327.       ddb = b[i][jprev] + b[i][jnext] + b[iprev][j] + b[inext][j] - 4 * bval;
  328.  
  329.       da[i][j] = ka * (16 - aval * bval) + Diff1 * dda;
  330.       db[i][j] = ka * (aval * bval - bval - c[i][j]) + Diff2 * ddb;
  331.     }
  332.   }
  333.  
  334.   /* affect change */
  335.  
  336.   for (i = 0; i < xsize; i++)
  337.     for (j = 0; j < ysize; j++) {
  338.       a[i][j] += (speed * da[i][j]);
  339.       b[i][j] += (speed * db[i][j]);
  340.       if (b[i][j] < 0)
  341.     b[i][j] = 0;
  342.     }
  343. }
  344.  
  345.  
  346. /******************************************************************************
  347. Calculate semi-stable equilibria.
  348. ******************************************************************************/
  349.  
  350. semi_equilibria()
  351. {
  352.   int i,j;
  353.   float ainit,binit;
  354.   float cinit,dinit,einit;
  355.  
  356.   ainit = binit = cinit = dinit = einit = 0;
  357.  
  358.   /* figure the values */
  359.  
  360.   switch (sim) {
  361.  
  362.     case STRIPES:
  363.       for (i = 0; i < xsize; i++) {
  364.  
  365.     ainit = p2 / (2 * p1);
  366.     binit = ainit;
  367.     cinit = 0.02 * ainit * ainit * ainit / p2;
  368.     dinit = ainit;
  369.     einit = ainit;
  370.     
  371.     for (j = 0; j < ysize; j++) {
  372.       a[i][j] = ainit;
  373.       b[i][j] = binit;
  374.       c[i][j] = cinit;
  375.       d[i][j] = dinit;
  376.       e[i][j] = einit;
  377.       ai[i][j] = 1 + frand (-0.5 * arand, 0.5 * arand);
  378.     }
  379.       }
  380.       break;
  381.  
  382.     case SPOTS:
  383.       for (i = 0; i < xsize; i++)
  384.     for (j = 0; j < ysize; j++) {
  385.       a[i][j] = a_steady;
  386.       b[i][j] = b_steady;
  387.       c[i][j] = beta_init + frand (-beta_rand, beta_rand);
  388.     }
  389.       break;
  390.  
  391.     default:
  392.       printf ("bad case in semi_equilibria\n");
  393.       break;
  394.   }
  395. }
  396.  
  397.  
  398. /******************************************************************************
  399. Switch for picking array to rescale.
  400. ******************************************************************************/
  401.  
  402. do_rescale (index,min,max)
  403.   int index;
  404.   float min,max;
  405. {
  406.   switch (index) {
  407.     case 1:
  408.       rescale_values (a, min, max);
  409.       break;
  410.     case 2:
  411.       rescale_values (b, min, max);
  412.       break;
  413.     case 3:
  414.       rescale_values (c, min, max);
  415.       break;
  416.     case 4:
  417.       rescale_values (d, min, max);
  418.       break;
  419.     case 5:
  420.       rescale_values (e, min, max);
  421.       break;
  422.     default:
  423.       printf ("bad switch in do_rescale: %d\n", index);
  424.       break;
  425.   }
  426. }
  427.  
  428.  
  429. /******************************************************************************
  430. Rescale values in array.
  431.  
  432. Entry:
  433.   values    - array to rescale
  434.   min_final - minimum value to map to
  435.   max_final - maximum value to map to
  436. ******************************************************************************/
  437.  
  438. rescale_values (values,min_final,max_final)
  439.   float values[MAX][MAX];
  440.   float min_final;
  441.   float max_final;
  442. {
  443.   int i,j;
  444.   float val;
  445.   float min =  1e20;
  446.   float max = -1e20;
  447.  
  448.   /* find minimum and maximum values */
  449.  
  450.   for (i = 0; i < xsize; i++)
  451.     for (j = 0; j < ysize; j++) {
  452.       if (values[i][j] < min)
  453.     min = values[i][j];
  454.       if (values[i][j] > max)
  455.     max = values[i][j];
  456.     }
  457.  
  458.   if (min == max) {
  459.     min = max - .001;
  460.     max = min + .002;
  461.   }
  462.  
  463.   /* rescale the values */
  464.  
  465.   for (i = 0; i < xsize; i++)
  466.     for (j = 0; j < ysize; j++) {
  467.       val = (values[i][j] - min) / (max - min);
  468.       val = min_final + val * (max_final - min_final);
  469.       values[i][j] = val;
  470.     }
  471. }
  472.  
  473.  
  474. /******************************************************************************
  475. Display the activator.
  476. ******************************************************************************/
  477.  
  478. show(values)
  479.   float values[MAX][MAX];
  480. {
  481.   int i,j;
  482.   float output;
  483.   float min =  1e20;
  484.   float max = -1e20;
  485.  
  486.   /* find minimum and maximum values */
  487.  
  488.   for (i = 0; i < xsize; i++)
  489.     for (j = 0; j < ysize; j++) {
  490.       if (values[i][j] < min)
  491.     min = values[i][j];
  492.       if (values[i][j] > max)
  493.     max = values[i][j];
  494.     }
  495.  
  496.   if (min == max) {
  497.     min = max - 1;
  498.     max = min + 2;
  499.   }
  500.  
  501.   printf ("min max diff: %f %f %f\n", min, max, max - min);
  502.  
  503.   /* display the values */
  504.  
  505.   for (i = 0; i < xsize; i++)
  506.     for (j = 0; j < ysize; j++) {
  507.       output = (values[i][j] - min) / (max - min);
  508.       output = output * 255.0;
  509.       writepixel (i, j, (int) output);
  510.     }
  511.  
  512.   flushbuffers();
  513. }
  514.  
  515.  
  516. /******************************************************************************
  517. Pick a random number between min and max.
  518. ******************************************************************************/
  519.  
  520. float frand(min, max)
  521.   float min,max;
  522. {
  523.   return (min + drand48() * (max - min));
  524. }
  525.  
  526.