// simple, 1-box NPZ and NPZD models, no light dependence, no sinking, no physics.
// included as examples with ASTroCAT.

class SimpleNPZ extends Ecosystem {

  Statevar N,P,Z,totalN; 
  Flux Fuptake, Fingest, Fpmort, Fzmort, Fexcret;
  FreeParam mu0, I0, ks, Ks, m, mz, epsilon;
  
  // system definition ----------------------
  
  SimpleNPZ() {
    N = addVar("nutrients","N","µM");
    P = addVar("phytoplankton","P","µM");
    P.sinks = true;
    Z = addVar("zooplankton","Z","µM");
    totalN = addGhostVar("total nitrogen","T","µM");
    
    Fuptake = addFlux("uptake","uptake", N, P);
    Fingest = addFlux("ingestion","ingest", P, Z);
    Fpmort = addFlux("P mortality","Pmort", P, N);
    Fzmort = addFlux("Z mortality","Zmort", Z, N);
    Fexcret = addFlux("excretion","excret", Z, N);

    mu0 = addParam("P max growth","mu0","1/d",2,0,2);
    ks = addParam("nutrient half-saturation","ks","µM",1, 0.1, 5);
    m = addParam("P mortality","m","1/d",0.1, 0, 0.4);
    I0 = addParam("Z max growth","I0","1/d", 3, 0, 5);
    Ks = addParam("grazing half-saturation","Ks","µM",3, 0, 4);
    epsilon = addParam("Z growth efficiency","epsilon","",0.3, 0, 1);
    mz = addParam("Z mortality","mz","1/d",0.4, 0, 1.0);
  }
  
  // grazing and mortality are linear: simple, prone to oscillations
  void calcFluxes() {
    for (int i=0; i<Ncells; i++) {
      Fuptake.current[i] = mu0.current * N.current[i] / (ks.current + N.current[i]) * P.current[i];
      Fingest.current[i] = I0.current * P.current[i] / (Ks.current + P.current[i]) * Z.current[i];
      Fpmort.current[i] = m.current * P.current[i];
      Fzmort.current[i] = mz.current * Z.current[i];
      Fexcret.current[i] = (1 - epsilon.current) * Fingest.current[i];
      totalN.current[i] = N.current[i] + P.current[i] + Z.current[i];
    }
  }
  
}


class SimpleNPZ_NetworkView extends NetworkView {
  
  SimpleNPZ eco;
  
  SimpleNPZ_NetworkView(Ecosystem eco) {
    this.eco = (SimpleNPZ)eco;
  }

  void define() {
    varScale = 50;
    fluxScale = 4;
    float theta0 = -120;
    addDot_polar(eco.Z, 0.7, theta0-240, styles.medRedColor);
    addDot_polar(eco.P, 0.7, theta0-120, styles.ltGreenColor);
    addDot_polar(eco.N, 0.7, theta0, styles.dkBlueColor);
    
    addArrow(eco.Fuptake, styles.ltGreenColor, 120);
    addArrow(eco.Fingest, styles.medRedColor, 120);
    addArrow(eco.Fexcret, styles.medBlueColor, 120);
    addArrow(eco.Fpmort, styles.medBlueColor, 120);
    addArrow(eco.Fzmort, styles.medBlueColor, 120);
  }
      
}




// ---------------------------------------------------------------------------------------------------------




class SimpleNPZD extends Ecosystem {
  
  Statevar N,P,Z,D,totalN; 
  Flux Fuptake, Fingest, Fegest, Fpmort, Fzmort, Fexcret, Fremin, Fsupply, Fsink;
  FreeParam S, mu0, ks, alpha, attSW, attP, m, I0, Ks, epsilon, fegest, zeta, Rremin, wsinkD;
  float nominalBoxDepth = 20; // used in converting the D sinking velocity to a fraction lost per day
  
  // system definition ----------------------
  
  SimpleNPZD() {
    N = addVar("nutrients","N","µM");
    P = addVar("phytoplankton","P","µM");
    Z = addVar("zooplankton","Z","µM");
    D = addVar("detritus","D","µM");
    D.sinks = true;
    totalN = addGhostVar("total nitrogen","T","µM");

    Fsupply = addFlux("N supply","supply", null, N);
    Fuptake = addFlux("uptake","uptake", N, P);
    Fingest = addFlux("ingestion","ingest", P, Z);
    Fegest = addFlux("egestion","egest", P, D);
    Fexcret = addFlux("excretion","excret", Z, N);
    Fpmort = addFlux("P mortality","Pmort", P, D);
    Fzmort = addFlux("Z mortality","Zmort", Z, D);
    Fremin = addFlux("remineralization","remin", D, N);
    Fsink = addFlux("D sinking","Dsink",D,null);

    S = addParam("N supply rate","S","µM/d",1,0,2);
    S.highlight = true;
    mu0 = addParam("P max growth","mu0","1/d",2,0,3);
    ks = addParam("nutrient half-saturation","ks","µM",1, 0.1, 5);
    m = addParam("P mortality","m","1/d",0.1, 0, 0.4);
    I0 = addParam("Z max growth","I0","1/d", 3, 0, 6);
    Ks = addParam("grazing half-saturation","Ks","µM",3, 0.1, 10);
    Ks.logScale = true;
    epsilon = addParam("Z growth efficiency","epsilon","", 0.3, 0, 1);
    fegest = addParam("egestion fraction","fegest","",0.5, 0, 1);
    zeta = addParam("Z quadratic mortality","zeta","1/µM/d",2.0, 0, 4.0);
    Rremin = addParam("remineralization rate","remin","1/d",0.1, 0.01, 1);
    Rremin.logScale = true;
    wsinkD = addParam("D sinking rate","wsinkD","m/day",10,0,100);
  }
  
  void calcFluxes() {
    for (int i=0; i<Ncells; i++) {
      Fsupply.current[i] = S.current;
      float nutlim = N.current[i] / (ks.current + N.current[i]);
      Fuptake.current[i] = mu0.current * nutlim * P.current[i];
      float Fgraz = I0.current * sq(P.current[i]) / (sq(Ks.current) + sq(P.current[i])) * Z.current[i];
      Fingest.current[i] = epsilon.current * Fgraz;
      Fegest.current[i] = (1 - epsilon.current) * fegest.current * Fgraz;
      Fexcret.current[i] = (1 - epsilon.current) * (1 - fegest.current) * Fgraz;
      Fingest.current[i] += Fexcret.current[i]; // because excretion is written as a flux from Z to N
      Fpmort.current[i] = m.current * P.current[i];
      Fzmort.current[i] = zeta.current * sq(Z.current[i]);
      Fremin.current[i] = Rremin.current * D.current[i];
      Fsink.current[i] = wsinkD.current / nominalBoxDepth * D.current[i]; // pretty rough representation of sinking
      totalN.current[i] = N.current[i] + P.current[i] + Z.current[i] + D.current[i];      
    }
  }

}

class SimpleNPZD_NetworkView extends NetworkView {
  
  SimpleNPZD eco;
  
  SimpleNPZD_NetworkView(Ecosystem eco) {
    this.eco = (SimpleNPZD)eco;
  }

  void define() {
    varScale = 50;
    fluxScale = 4;
    float theta0 = -120;
    Dot Ddot = addDot_polar(eco.D, 0.7, theta0-270, styles.ltBlueColor);
    addDot_polar(eco.Z, 0.7, theta0-180, styles.medRedColor);
    addDot_polar(eco.P, 0.7, theta0-90, styles.ltGreenColor);
    Dot Ndot = addDot_polar(eco.N, 0.7, theta0, styles.dkBlueColor);
    
    addArrow(eco.Fsupply, Ndot.xcen, Ndot.ycen-0.3, styles.dkBlueColor, 0);
    addArrow(eco.Fuptake, styles.ltGreenColor, 90);
    addArrow(eco.Fingest, styles.medRedColor, 90);
    addArrow(eco.Fegest, styles.medBlueColor, 30);
    addArrow(eco.Fexcret, styles.medBlueColor, 30);
    addArrow(eco.Fpmort, styles.medBlueColor, -30);
    addArrow(eco.Fzmort, styles.medBlueColor, 90);
    addArrow(eco.Fremin, styles.dkBlueColor, 90);
    addArrow(eco.Fsink, Ddot.xcen, Ddot.ycen-0.3, styles.medBlueColor, 0);
  }
  
}