Model of plankton dynamics

makefile

comp=clang
opt= -O3 -DWall -lgsl -lgslcblas -DHAVE_INLINE -lm
src=src.c
out=lagoon

$(out): $(src)
	$(comp) $(src) -o $(out) $(opt)

simuls: $(out)
	./$(out)

clean:
	rm output/*.txt

src.c

// Compile with make
// Change the comp variable to gcc or clang

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf_log.h>
#include <gsl/gsl_sf_exp.h>

#define FNSIZE 128

double Pow10(double n)
{
   double res;
   res = gsl_sf_exp(n * gsl_sf_log(10));
   return res;
}

int main(int argc, char *argv[])
{
   // RNG init
   gsl_rng *rng = gsl_rng_alloc(gsl_rng_taus2);
   gsl_rng_set(rng, time(NULL));

   FILE *sim_out;

   const double SCALING = 0.1;

   const double substeps = 100.0;
   double current_substep;
   
   const int timesteps = 30;
   int current_timestep;

   const int replicates = 20;
   int current_replicate;

   const double a = 0.12;
   const double b = 0.12;

   const double qp = 0.06;
   const double qv = 0.05;

   const double rv = 1.7;
   const double rp = 0.4 * rv;
   double drv = 0.0;

   const double init_density = 5.0;

   double V, P;
   double dV, dP;

   double LNmu;
   for(LNmu = 0.10 ; LNmu <= 1.1 ; LNmu += 0.5)
   {
      //double rLNmu = Pow10(LNmu);
      double rLNmu = LNmu;
      for(current_replicate = 0 ; current_replicate < replicates ; current_replicate++)
      {
         // Output file
         char tfname[FNSIZE];
         snprintf(tfname, sizeof(char) * FNSIZE, "output/h%.6f-r%d.txt", rLNmu, current_replicate);
         sim_out = fopen(tfname, "w");
         // Start simul
         V = init_density;
         P = init_density;
         for(current_timestep = 0 ; current_timestep < timesteps ; current_timestep++)
         {
            // Draw drv for this timestep
            drv = rv + gsl_ran_lognormal(rng, 0.0, rLNmu) * SCALING;
            for(current_substep = 0.0 ; current_substep < substeps ; current_substep += 1.0)
            {
               dV = (V * drv - qv * V * V -a * V * P) * (1.0 / substeps);
               dP = (P * rp - qp * P * P + a * b * V * P) * (1.0 / substeps);
               V += dV;
               P += dP;
               V = V > 0.0 ? V : 0.0;
               P = P > 0.0 ? P : 0.0;
            }
            if(current_timestep >= 5) fprintf(sim_out, "%d %.5f %.5f %.6f\n", current_timestep, V, P, drv);
         }
         fclose(sim_out);
      }
   }

   // Free and return
   gsl_rng_free(rng);
   return EXIT_SUCCESS;
}