class ASTroCAT extends DiverseEcosystem {
  
  // the ASTroCAT ecosystem model itself.
  //
  // imposed N supply, spectrum of P, spectrum of Z, and that's it. No light dependence.
  
  Statevar Ptot, Ztot, N;
  Statevar[] P, Z;
  
  Flux[][] FgrazP; // one for each P-Z pair
  Flux[] Fegest; // summed over Z for each P
  Flux[] Fexcret; // summed over P for each Z
  Flux[] Fuptake, Fpmort, Fzmort;
  Flux Fsupply;
  
  FreeParam epsilon, fegest, supply; // system params
  FreeParam mu0, m, I0, zeta;  // rate params
  FreeParam ks, K0; // shape params
  FreeParam xpreyopt, dxprey; // grazing pref params
  float[][] phiP; // grazing preferences P -> Z
  boolean doAnnualCycle = false;
  
  int NP, NZ;
  
  // system definition ----------------------
  
  
  
  ASTroCAT() {construct(10);}
  
  ASTroCAT(int NP) {construct(NP);}
    
  void construct(int NP) {
    this.NP = NP;
    this.NZ = NP;
    
    float PMinEsd = 1; // size ranges in microns
    float PMaxEsd = 20;
    
    // state variables
    N = addVar("nutrients","N","µM");
    P = addDiverseVarSet("phytoplankton","P","µM",NP,PMinEsd, PMaxEsd);
    Ptot = addGhostVar("total P","Ptot","µM");
    Z = addDiverseVarSet("zooplankton","Z","µM",NZ,1,2);
    for (int i=0; i<NP; i++) ((DiverseVar)Z[i]).esd = 2.16 * pow(((DiverseVar)P[i]).esd, 1.79);
    Ztot = addGhostVar("total Z","Ztot","µM");
    
    // minimum biomass values: always protect a little bit of seed stock
    for (int i=0; i<NP; i++) P[i].min = 0.01;
    for (int j=0; j<NZ; j++) Z[j].min = 0.01;
    
    // fluxes
    Fsupply = addFlux("supply","supply",null,N);
    Fuptake = addFluxSet("uptake","uptake", N, P);
    FgrazP = addFluxSet2D("grazing on P","grazP",P, Z);
    Fegest = addFluxSet("egestion","egest", Z, null);
    Fexcret = addFluxSet("excretion","excret", Z, N);
    Fpmort = addFluxSet("P mortality","Pmort", P, null);
    Fzmort = addFluxSet("Z mortality","Zmort", Z, null);

    // global params
    supply = addParam("N supply rate", "supply","µM/d", 1, 0.01, 4);
    supply.logScale = true;
    supply.highlight = true;
    FreeParam p = addParam("annual cycle on/off","annCycle","", 0, 0, 1);
    p.quantize(1);
    p.highlight = true;    
    epsilon = addParam("Z growth efficiency","epsilon","", 0.33, 0, 1);
    fegest = addParam("egestion fraction","fegest","",0.33, 0, 1); // a fraction of total grazing, not of the losses
    m = addParam("P mortality as fraction of mu0","m","", 0.1, 0, 1);
    dxprey = addParam("prey size preference width","dxprey","", 0.25, 0.01, 1);
    dxprey.logScale = true;
    
    // diverse params: defined using a scale size of 1 µm for easy reading
    // std dev within a size class is defined based on the sources cited in the astrocat paper, but not used in astrocat 1.0
    //   rate params
    mu0 = addDiverseParam("P max growth","mu0","1/d",                   2.6,     0,   3,   3.8,    1,  -0.45); // initial      min      max     std dev (*,/) within a size class      allo scale size     exponent for allometry
    I0 = addDiverseParam( "Z max grazing","I0","1/d",                    26,     0,  30,   1.9,    1,  -0.4);
    zeta = addDiverseParam("Z quadratic mortality","zeta","1/µM/d",       1,     0,  10,     1,    1,   0);
    zeta.highlight = true;
    //  functional-response shape params
    ks = addDiverseParam("nutrient half-saturation","ks","µM",          0.1,   0.1,   5,     2,    1,   1);
    K0 = addDiverseParam("grazing half-saturation","K0","µM",             3,   0.3,  30,     3,    1,   0);
    K0.logScale = true;
    //  grazing preference params
    xpreyopt = addDiverseParam("optimal prey size","xpreyopt","",      0.65,   0.5,  30,     1,    1,   0.56); // has allometry but no stochasticity; noise in prey prefs is handled separately

    phiP = new float[NP][NZ]; // prey preferences: internal array
    addParam("max preference","maxPref","", 10, 0, 10);

    p = addParam("prefs allometry on/off","prefsAllo","",1, 0, 1);
    p.quantize(1);
    p.highlight = true;
    p = addParam("other allometry on/off","otherAllo","",1, 0, 1);
    p.quantize(1);
    p.highlight = true;    
    p = addParam("prefs stochasticity","prefsRand","",0, 0, 1); // other stochasticity (in mu0, I0, for example) is omitted in this version of the model
    p.quantize(1);
    p.highlight = true;        
 
   // set a better scale size, so the monotype version doesn't run with 1 µm cells
    float xP = sqrt(PMinEsd*PMaxEsd); // geometric means of the size ranges
    float xZ = 2.16 * pow(xP, 1.79); // set Z sizes such that each P is optimal prey for one Z
    ((DiverseParam)mu0).initial *= pow(xP / ((DiverseParam)mu0).alloScaleSize, ((DiverseParam)mu0).alloExp);
    ((DiverseParam)mu0).alloScaleSize = xP;
    ((DiverseParam)I0).initial *= pow(xZ / ((DiverseParam)I0).alloScaleSize,  ((DiverseParam)I0).alloExp);
    ((DiverseParam)I0).alloScaleSize = xZ;
    ((DiverseParam)xpreyopt).initial *= pow(xZ / ((DiverseParam)xpreyopt).alloScaleSize,  ((DiverseParam)xpreyopt).alloExp);
    ((DiverseParam)xpreyopt).alloScaleSize = xZ;
    ((DiverseParam)ks).initial *= pow(xP / ((DiverseParam)ks).alloScaleSize,  ((DiverseParam)ks).alloExp);
    ((DiverseParam)ks).alloScaleSize = xP;
    // adjust the names of the allometric params so the user can see what the sliders control
    mu0.name += " at " + round(xP) + " µM";
    I0.name += " at " + round(xZ) + " µM";
    ks.name += " at " + round(xP) + " µM";
    xpreyopt.name += " at " + round(xZ) + " µM";
  }
  
  
  
  void calcFluxes() {
    for (int k=0; k<Ncells; k++) { // for each spatial cell...
      // supply
      Fsupply.current[k] = supply.current;
      if (doAnnualCycle) Fsupply.current[k] *= 0.5*(1+cos(time.current[0] * TWO_PI/365.));
      // phyto growth
      for (int i=0; i<NP; i++) {
        float mu0i = getCurrent(mu0, P[i]);
        float ksi = getCurrent(ks, P[i]);
        float nutlim = N.current[k] / (ksi + N.current[k]);
        Fuptake[i].current[k] = mu0i * nutlim * P[i].current[k];
      }
      
      // grazing
      //     assemble all diet elements
      float[][] PscaledAsFood = new float[NP][NZ];
      float[] totalFood = new float[NZ];
      for (int j=0; j<NZ; j++) {
        float K0j = getCurrent(K0, Z[j]);
        totalFood[j] = 0;
        for (int i=0; i<NP; i++) {
          PscaledAsFood[i][j] = phiP[i][j] / K0j * P[i].current[k];
          totalFood[j] += PscaledAsFood[i][j];
        }
      }
      //     initialize excretion and egestion fields
      for (int i=0; i<NP; i++) Fegest[i].current[k] = 0;
      for (int j=0; j<NP; j++) Fexcret[j].current[k] = 0;
      //     calculate fluxes
      for (int j=0; j<NZ; j++) {
        float I0j = getCurrent(I0, Z[j]);
        for (int i=0; i<NP; i++) {
          FgrazP[i][j].current[k] = I0j * Z[j].current[k] * PscaledAsFood[i][j] / (1. + totalFood[j]);
          Fegest[j].current[k] += fegest.current * FgrazP[i][j].current[k];
          Fexcret[j].current[k] += (1.0 - fegest.current - epsilon.current) * FgrazP[i][j].current[k];          
        }
      }
      
      // mortality
      for (int i=0; i<NP; i++) {
        float mi = getCurrent(m, P[i]);
        float mu0i = getCurrent(mu0, P[i]);
        Fpmort[i].current[k] = mi * mu0i * P[i].current[k];
      }
      for (int j=0; j<NZ; j++) {
        float zetaj = getCurrent(zeta, Z[j]);
        Fzmort[j].current[k] = zetaj * Z[j].current[k] * Ztot.current[k];
      }
      
      // diagnostics
      Ptot.current[k] = 0;
      for (int i=0; i<NP; i++) Ptot.current[k] += P[i].current[k];
      Ztot.current[k] = 0;
      for (int j=0; j<NZ; j++) Ztot.current[k] += Z[j].current[k];
    }
    
  }



  void updatePrivateParams() {
    doAnnualCycle = (getParamVal("annCycle")==1);
    
    // turn off stochasticity in rates and half-saturations
    ((DiverseParam)mu0).randMultiplier = 0;
    ((DiverseParam)ks).randMultiplier = 0;
    ((DiverseParam)I0).randMultiplier = 0;
    ((DiverseParam)K0).randMultiplier = 0;
    ((DiverseParam)zeta).randMultiplier = 0;
    
    // turn allometry on and off, in response to the user's choice
    float r = getParamVal("otherAllo");
    ((DiverseParam)mu0).doAllo = (r==1);
    ((DiverseParam)ks).doAllo = (r==1);
    ((DiverseParam)I0).doAllo = (r==1);
    ((DiverseParam)K0).doAllo = (r==1);
    ((DiverseParam)zeta).doAllo = (r==1);
    
    // set grazing preferences, including allometry and stochasticity
    r = getParamVal("prefsRand");
    boolean doAllo = (getParamVal("prefsAllo")==1);
    float maxPref = getParam("maxPref").current;
    for (int j=0; j<NZ; j++) {
      float xpreyoptj = getCurrent(xpreyopt, Z[j]);
      float dxpreyj = getCurrent(dxprey, Z[j]);
      for (int i=0; i<NP; i++) {
        float xpreyi = ((DiverseVar)P[i]).esd;
        float R = randn(((DiverseVar)P[i]).randSignature,((DiverseVar)Z[j]).randSignature);
        phiP[i][j] = exp(r*R);
        if (doAllo) phiP[i][j] *= exp(-sq((log10(xpreyi)-log10(xpreyoptj)) / dxpreyj));
        phiP[i][j] = constrain(phiP[i][j], 0, maxPref);
      }
    }
    
  }

}



class ASTroCAT_NetworkView extends DiverseNetworkView {
  
  ASTroCAT eco;
  
  ASTroCAT_NetworkView(Ecosystem eco) {
    this.eco = (ASTroCAT)eco;
  }
  
  void define() {
    varScale = 10 * (40.0/eco.NP);
    fluxScale = 0.5 * (40.0/eco.NP);

    float R = 0.8;
    float xcen = 0.2;
    float ycen = 0.1;
    int margin = 5;

    float[] px = new float[eco.NP];
    float[] py = new float[eco.NP];
    float[] zx = new float[eco.NZ];
    float[] zy = new float[eco.NZ];
      
    float[][] vertices = fluxArc(xcen-R, ycen-R, xcen-R+2*R*cos(PI/8), ycen-R+2*R*sin(PI/8), 30, eco.NZ+margin);
    for (int j=0; j<eco.NZ; j++) {
      zx[j] = vertices[j+margin][0];
      zy[j] = vertices[j+margin][1];
      Dot d = addDot(eco.Z[j], zx[j], zy[j], transparency(styles.dkBrownColor,0.8));
      d.name = nicenum(round(((DiverseVar)eco.Z[j]).esd));
    }
    vertices = fluxArc(xcen-R, ycen-R, xcen-R+2*R*cos(HALF_PI-PI/8), ycen-R+2*R*sin(HALF_PI-PI/8), 30, eco.NP+margin);
    for (int i=0; i<eco.NP; i++) {
      px[i] = vertices[i+margin][0];
      py[i] = vertices[i+margin][1];
      Dot d = addDot(eco.P[i], px[i], py[i], transparency(styles.medGreenColor,0.8));
      d.name = nicenum(round(((DiverseVar)eco.P[i]).esd));
    }
    addDot(eco.N, xcen-R, ycen-R, styles.dkBlueColor);
    
    // supply arrow
    Arrow a = addArrow(eco.Fsupply, xcen-2*R, ycen-R, styles.dkBlueColor, 0);
    a.arrowheadPos = 0.8;

    // regeneration and loss arrows
    addArrowSet(eco.Fexcret, styles.ltGrayColor, 150);
    for (int j=0; j<eco.NZ; j++) {
      a = addArrow(eco.Fzmort[j], zx[j]+0.4, zy[j], styles.medGrayColor, 0);
      a.setFrom(zx[j], zy[j]);
    }
    for (int i=0; i<eco.NP; i++) {
      a = addArrow(eco.Fpmort[i], px[i]+0.4, py[i], styles.medGrayColor, 0);
      a.setFrom(px[i], py[i]);
    }
    
    // growth & grazing arrows
    Arrow[] aa = addArrowSet(eco.Fuptake, styles.ltGreenColor, 120);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
    aa = addArrowSet(eco.FgrazP, styles.dkBrownColor, 90);
    for (int i=0; i<aa.length; i++) aa[i].arrowheadPos = 1;
  }
  
  void updateArrows() {
    for (int i=0; i<arrows.length; i++) {
      arrows[i].update(fluxes[i].current[currentCell]);
    }
    if (arrows.length > 400) {
      arrows[0].draw(); // draw the supply arrow only if there's a lot of arrows
    } else {
      for (int i=0; i<arrows.length; i++) arrows[i].draw();
    }
  }
  
}