class dynamicalSystem { // -----------------------------------------------------------
  String name,shortname;
  String timeUnits,varUnits;
  int length, nvars, nfluxes, nparams;
  flux[] fluxes;
  statevar[] vars;
  freeParam[] params;

  float[][] varhistory,fluxhistory,paramhistory;
  float[] historytime;
  int nhistend = -1;
  boolean wrapHistory = false;
  float tcurrent = 0;
  float dtlast = nan;
  
  float modelTimePerSecond = 1;
  float tolerance = 1e-4;
  boolean adjustFluxesAtMaxMin = false; // true = adjust fluxes so statevars don't exceed their max/mins;
                                        // false = constrain statevars but pretend full fluxes were applied
  
  float dot_stdsize = 0.2; // relative to diagram width
  float dot_minsize = 0.01;
  float dot_rigidity = 1; // how linear the scaling is:  >> 1, hard to change dot size; << 1, close to linear
  float line_stdwidth = 4;
  float line_minwidth = 0.25;
  float line_maxwidth = 10;
  float line_rigidity = 1;
  int arcResolution = 40;
  boolean showFluxes = true;
  boolean showDetails = false;
  boolean showDetailsButton = false;
  boolean showVarnames = false;
  boolean showHistory = true;
  boolean showNewtork = true;
  boolean showVarDiagnostics = false;

  int[] networkRect = new int[4];
  int[] historyRect = new int[4];
  
 
 
  void allocate(int N, int maxvars, int maxfluxes, int maxparams) {
    length = N;
    fluxes = new flux[maxfluxes];
    nfluxes = 0;
    vars = new statevar[maxvars];
    nvars = 0;
    params = new freeParam[maxparams];
    nparams = 0;
    varhistory = new float[maxvars][N];
    //fluxhistory = new float[maxfluxes][N]; // fluxhistory commented out here and in savetohistory(), so that the amount of
                                               // history storage goes like (nvars), not (nvars^2)
    paramhistory = new float[maxparams][N];
    historytime = new float[N];
  }
  
  dynamicalSystem clone() {
    println("warning: no clone() method has been defined for this system");
    return null;
  }
  
  
  int addFlux(String thename, String theshortname, int thefrom, int theto, color thecolor, float curv) {
    flux newflux;
    newflux = new flux(length,nan);
    newflux.define(thename, theshortname, thefrom, theto, thecolor, curv);
    nfluxes++;
    fluxes[nfluxes-1] = newflux;
    return nfluxes-1;
  }
  
  int addVar(String thename, String theshortname, float initialval, color thecolor, float theX, float theY) {
    statevar newvar;
    newvar = new statevar(length,nan);
    newvar.define(thename, theshortname, initialval, thecolor, theX, theY);
    nvars++;
    vars[nvars-1] = newvar;
    return nvars-1;
  }
  
  int addParam(String thename, String theshortname, float initialval, float minval, float maxval, color thecolor) {
    freeParam newparam;
    newparam = new freeParam(length,nan);
    newparam.define(thename, theshortname, initialval, minval, maxval, thecolor);
    nparams++;
    params[nparams-1] = newparam;
    return nparams-1;
  }
  
  int addEndpoint(String thename, String theshortname, float theX, float theY) { // special kind of statevar
    statevar newvar;
    newvar = new statevar(length,nan);
    newvar.define(thename, theshortname, nan, nan, theX, theY);
    nvars++;
    vars[nvars-1] = newvar;
    vars[nvars-1].constant = true;
    vars[nvars-1].display = false;
    return nvars-1;
  }
  
  int addConstant(String theshortname, float val) {
    return addConstant(theshortname,theshortname,val);
  }
  int addConstant(String thename, String theshortname, float initialval, float minval, float maxval, color thecolor) {
    // if you addConstant with the arguments you'd use for addParam, the extra ones are ignored
    return addConstant(thename, theshortname, initialval);
  }
  int addConstant(String thename, String theshortname, float val) { // special kind of param
    int id = addParam(thename,theshortname,val,val,val,color(0));
    params[id].frozen = true;
    return id;
  }
  
  
  int find(String sn) {
  // returns the index into fluxes, vars, or params that matches the shortname "sn".
  // warning: the first match is returned; there's no checking for duplicates. So make sure there
  // aren't any repetitions of shortnames among the 3 arrays.
    int a = -1;
    int i = -1;
    while ((a==-1) && (i<nfluxes-1)) {i++;  if (sn==fluxes[i].shortname) {a=i;}}
    i = -1;
    while ((a==-1) && (i<nvars-1)) {i++;  if (sn==vars[i].shortname) {a=i;}}
    i = -1;
    while ((a==-1) && (i<nparams-1)) {i++;  if (sn==params[i].shortname) {a=i;}}
    return a;
  }
  
  
  float readout(String sn) {
  // returns the value of the var, flux, or param that matches the shortname "sn".
    float V = nan;
    int i = -1;
    while ((V==nan) && (i<nfluxes-1)) {i++; if (sn==fluxes[i].shortname) {V=fluxes[i].current;}}
    i = -1;
    while ((V==nan) && (i<nvars-1)) {i++; if (sn==vars[i].shortname) {V=vars[i].current;}}
    i = -1;
    while ((V==nan) && (i<nparams-1)) {i++; if (sn==params[i].shortname) {V=params[i].current;}}
    return V;
  }
  
  
  void wakeup() {
    reset_varsfluxes();
    reset_params();
  }
  
  void select() {}
  
  
  float[] addsomeflux(float[] V, float[] F, float dt) {
  // performs the basic operation Vout = V + F*dt. This is what calc() did in previous versions, but now
  // calc() is a wrapper that handles the adaptive stepsize part.
    int jto,jfrom; 
    float dV;
    float[] Vout = new float[nvars];
    for (int j=0; j<nvars; j++) {Vout[j] = V[j];}
    for (int i=0; i<nfluxes; i++) {
      dV = F[i] * dt;
      jto = fluxes[i].to;
      jfrom = fluxes[i].from;

      if ((adjustFluxesAtMaxMin) && ((!vars[jto].constant) && (!vars[jfrom].constant))) {
        // adjust flux so that vars don't exceed the set max & min: necessary for mass conservation (have to tell _to_ if _from_ is limiting the flux and vice versa)
        if (dV > 0) {
          dV = min(dV, vars[jto].max - Vout[jto]);
          dV = min(dV, Vout[jfrom] - vars[jfrom].min);
        } else {
          dV = max(dV, vars[jto].min - Vout[jto]);
          dV = max(dV, Vout[jfrom] - vars[jfrom].max);
        }
        dV = constrain(dV, fluxes[i].min * dt, fluxes[i].max * dt);
      }

      if (!vars[jto].constant) {
        Vout[jto] += dV;
        Vout[jto] = constrain(Vout[jto],  vars[jto].min, vars[jto].max);
      }
      if (!vars[jfrom].constant) {
        Vout[jfrom] -= dV;
        Vout[jfrom] = constrain(Vout[jfrom],  vars[jfrom].min, vars[jfrom].max);        
      }     
    }  
    return Vout;
  }
  
  
  // these are needed to make romsBioSlice work. Not elegant.
  void setInitialConditions() {}
  void setInitialConditions(float someInputVal) {}
  float totalP() {return 0;}
  float totalNitrogen() {return 0;}
  

  void adjustments_after_calc_step() { // placeholder for things like recalculating total nitrogen and statevar scales 
  }
  
  float[] CALCFLUXES(float[] vslice) { // placeholder for the routine where the system-specific ODEs go
    return null;
  }
  
  int calc(float DT) {   // integrates the system from tcurrent to tcurrent+DT.              
    float[] Vcurrent,Pcurrent,fslice;
    fslice = new float[nfluxes]; Vcurrent = new float[nvars]; Pcurrent = new float[nparams];
    float[] Vfull,Vhalf;
    Vfull = new float[nvars]; Vhalf = new float[nvars];
    float err,thiserr,DTsofar,dt;
    int errorMsg = 0;
        
    for (int i=0; i<nvars; i++) {Vcurrent[i] = vars[i].current;}
    // DT is the total distance we have to integrate across; dt is the current stepsize along the way 
    DTsofar = 0;
    if (dtlast==nan) {dtlast = fluxes[0].scale*tolerance;}
    dt = dtlast;
    while (DTsofar < DT) {
      if (DT-DTsofar < dt) {dt = DT-DTsofar;} // don't integrate _too_ far
      // take a full step...
      fslice = CALCFLUXES(Vcurrent);
      Vfull = addsomeflux(Vcurrent,fslice,dt);
      // now take two half steps...
      Vhalf = addsomeflux(Vcurrent,fslice,dt/2);
      fslice = CALCFLUXES(Vhalf);
      Vhalf = addsomeflux(Vhalf,fslice,dt/2);
      // was this good enough?
      err = 0;
      int worstVar = 0;
      for (int i=0; i<nvars; i++) {
        thiserr = abs(Vhalf[i] - Vfull[i])/vars[i].scale;
        if (thiserr>err) {
          err = thiserr;
          worstVar = i;
        }
      }
      if ((err <= tolerance) || (dt < tolerance*DT)) { // successful
                             // (or dt is so small we'd better call it successful and move on)
        DTsofar = DTsofar + dt;
        Vcurrent = Vhalf; // note: haven't saved it to system.vars yet
        if (err>tolerance) {
          errorMsg = 1;
//          if (debug) {println(shortname + ": " + vars[worstVar].shortname + ": half full scale = " + Vhalf[worstVar] + "  " + Vfull[worstVar] + "  " + vars[worstVar].scale);}
        }
        dt = 1.2 * dt; // try a bigger step next time
      } else { // unsucessful
        dt = 0.5 * dt; // cut the stepsize in half and try again
      }
    }
    tcurrent = tcurrent + DT;  
    for (int i=0; i<nvars; i++) {vars[i].current = Vcurrent[i];} 
    fslice = CALCFLUXES(Vcurrent);
    for (int i=0; i<nfluxes; i++) {fluxes[i].current = fslice[i];}
    dtlast = dt; // for next time, a record of the last successful stepsize 
    
    adjustments_after_calc_step();
       
    return errorMsg;
    
  } // calc
  
  
  void savetohistory() {
    // at any given moment, the history record runs from index nhiststart to nhistend. If nhiststart comes after nhistend,
    // that means wrap around.
    // so add 1 to nhistend in a wraparound way, save the next frame there, and then figure out where nhiststart should be.
    
    nhistend++;
    if (nhistend==length) {
      nhistend = 0;
      wrapHistory = true;
    }
    for (int i=0; i<nvars; i++) {varhistory[i][nhistend] = vars[i].current;}
    //for (int i=0; i<nfluxes; i++) {fluxhistory[i][nhistend] = fluxes[i].current;}
    for (int i=0; i<nparams; i++) {paramhistory[i][nhistend] = params[i].current;}
    historytime[nhistend] = tcurrent;
  }
  
  
  void onmousepress() {
    for (int i=0; i<nparams; i++) {params[i].onMousePress();}
  }
  
  void update() {
    boolean needtorespond = false;
    setFont("med");
    for (int i=0; i<nparams; i++) { // update parameter controllers, if any
      params[i].update();
      if (params[i].current != params[i].previous) {needtorespond = true;}
    }
    if (needtorespond) {
      RESPOND_TO_PARAMCHANGE();
    }
    for (int i=0; i<nparams; i++) {
      params[i].previous = params[i].current;
    }
  }
  
  void showFluxArrows() {showFluxes = true;}
  void hideFluxArrows() {showFluxes = false;}
  void showGraphicDetails() {showDetails = true;}
  void hideGraphicDetails() {showDetails = false;}

  void RESPOND_TO_PARAMCHANGE() {} // placeholder for bookkeeping that occurs when parameters are changed by the user
  
  void reset_varsfluxes() {
    // resets system to initial state, but preserves params.
    for (int i=0; i<nfluxes; i++) {fluxes[i].current = fluxes[i].initial;}
    for (int i=0; i<nvars; i++) {vars[i].current = vars[i].initial;}
    tcurrent = 0;
    nhistend = -1;
    wrapHistory = false;
    savetohistory();
    dtlast = nan;
    adjustments_after_calc_step();
  }
  
  void reset_params() {
    for (int i=0; i<nparams; i++) {params[i].current = params[i].initial;}
    update();
  }

  // adjust dot & arrow sizes, positions, & visibility following an integration step, if the default method is insufficient
  void adjustments_before_drawing() {}
  
  
  void draw_network() {
    float w,r,dir,total;
    float xfrom,yfrom,xto,yto,xcen,ycen;
    float xscale,xoffset,yscale,yoffset;
    
    // conversion from relative coords -1..1 to screen coords
    xscale = networkRect[2]/2;
    xoffset = networkRect[0]+xscale;
    yscale = networkRect[3]/2;
    yoffset = networkRect[1]+yscale;

    // update dot sizes and arrow widths and set currentcolor to default color...
    for (int i=0; i<nvars; i++) {
      vars[i].radius = constrain(dot_stdsize * sqrt(vars[i].current/vars[i].scale), dot_minsize, 1);
      vars[i].currentcolor = vars[i].thecolor;
    }
    for (int i=0; i<nfluxes; i++) {
      fluxes[i].width = constrain(line_stdwidth * abs(fluxes[i].current/fluxes[i].scale), line_minwidth, line_maxwidth);
      fluxes[i].currentcolor = fluxes[i].thecolor;
      if (abs(fluxes[i].current) < tolerance*abs(fluxes[i].scale)) {fluxes[i].currentcolor = 255;}
    }
   
    // ...override default graphics settings if so desired
    adjustments_before_drawing();    
                    
    // draw flux arrows
    if (showFluxes) {
      for (int i=0; i<nfluxes; i++) {
        if (fluxes[i].current != nan) {
          xfrom = vars[fluxes[i].from].xc * xscale + xoffset;
          yfrom = -vars[fluxes[i].from].yc * yscale + yoffset;
          xto = vars[fluxes[i].to].xc * xscale + xoffset;
          yto = -vars[fluxes[i].to].yc * yscale + yoffset;
          drawonefluxarrow(i,xfrom,yfrom,xto,yto);
        }
      }
    }
    // draw statevar dots
    for (int i=0; i<nvars; i++) {
      if ((vars[i].display) && (vars[i].current != nan)) {
        r = vars[i].radius * xscale;
        xcen = vars[i].xc*xscale+xoffset;
        ycen = -vars[i].yc*yscale+yoffset;
        drawonevardot(i,xcen,ycen,r);
        if ((showVarDiagnostics) && (sqrt(sq(mouseX-xcen)+sq(mouseY-ycen))) <= r) { // write the value of the variable if the mouse is over it
          fill(0);
          noStroke();
          setFont("small");
          textAlign(RIGHT);
          text(""+round(vars[i].current*100)/100.0 + " " + varUnits, xcen-0.9*r, ycen+0.9*r);
        }
        if (showVarnames) {
          drawonevarlabel(i,xcen,ycen,r);      
        }  
      }
    } 
  }
  
  void drawonevardot(int i, float xcen, float ycen, float r) {
    noStroke();
    fill(vars[i].currentcolor);
    if (brightness(vars[i].currentcolor) > 240) {stroke(transparency(dkGrayColor,0.5)); strokeWeight(1);}
    ellipseMode(CENTER);
    ellipse(xcen,ycen,2*r,2*r);
  }
  
  void drawonevarlabel(int i, float xcen, float ycen, float r) {
    color col = vars[i].currentcolor;
    textAlign(CENTER);
    if (r>smallfontsize) {
      setFont("flexi");
      fill(lighten(lighten(lighten(col))));
      if (brightness(col) > 240) {fill(dkGrayColor);}
      textSize(r);
      textVAlign(vars[i].shortname,xcen,ycen,r,"middle");
    } else {
      setFont("small");
//      textVAlign(vars[i].shortname,xcen,ycen-r-1,smallfontsize,"bottom");
      float theta = atan2(vars[i].yc,vars[i].xc);
      textVAlign(vars[i].shortname,xcen+(vars[i].radius+smallfontsize)*cos(theta),ycen-(vars[i].radius+smallfontsize)*sin(theta),smallfontsize,"middle");
    }
  }
  
  void drawonefluxarrow(int i, float xfrom, float yfrom, float xto, float yto) {
    strokeWeight(fluxes[i].width);
    stroke(fluxes[i].currentcolor);
    float[][] vertices = flux_arc(xfrom,yfrom,xto,yto,fluxes[i].curvature,arcResolution);
    noStroke();
    fill(fluxes[i].thecolor);
    if (fluxes[i].current>0) {
      arrowhead(vertices,fluxes[i].arrowheadpos,fluxes[i].width*6,1);
    } else if (fluxes[i].current<0) {
      arrowhead(vertices,1-fluxes[i].arrowheadpos,fluxes[i].width*6,-1);
    }
  }
  
 
  void draw_history_setup(int n0, int n1) {
    float xscale,yscale,xoffset,yoffset;
    float timePrecision = 10;
    int n;
    
    fill(lighten(ltGrayColor));    noStroke();
    rect(historyRect[0], historyRect[1], historyRect[2], historyRect[3]);
    
    strokeWeight(1.5);    stroke(dkGrayColor);    noFill();
    beginShape();
    vertex(historyRect[0],historyRect[1]);
    vertex(historyRect[0],historyRect[1]+historyRect[3]);
    vertex(historyRect[0]+historyRect[2],historyRect[1]+historyRect[3]);
    endShape();
      
    int thetextsize = medfontsize;
    setFont("med");    fill(dkGrayColor);    noStroke();
    float t = (float)round(historytime[n0]*timePrecision)/timePrecision;
    textAlign(CENTER);
    textVAlign(str(t), historyRect[0], historyRect[1]+historyRect[3]+4, thetextsize, "top");
    t = (float)round(historytime[n1]*timePrecision)/timePrecision;
    textAlign(RIGHT);
    textVAlign(str(t), historyRect[0]+historyRect[2], historyRect[1]+historyRect[3]+4, thetextsize, "top");
    textAlign(LEFT);
    textVAlign(" " + timeUnits, historyRect[0]+historyRect[2], historyRect[1]+historyRect[3]+4, thetextsize, "top");
    textAlign(RIGHT);
    
    float total = 0;
    for (int i=0; i<nvars; i++) {if (!vars[i].constant) {total = total + vars[i].current;}}
    float totalrounded = (float)round(total*timePrecision)/timePrecision;
    textVAlign(str(totalrounded)+" "+varUnits+" ", historyRect[0],historyRect[1], thetextsize, "middle");
  }
  
  
  void draw_history_content(int n0, int n1) { // draw a history line for each parent dot
    for (int i=0; i<nvars; i++) { 
      if (vars[i].display) {
        float xscale = (float)historyRect[2]/constrain(historytime[n1]-historytime[n0],tolerance,1/tolerance);
        float xoffset = historyRect[0] - historytime[n0]*xscale;
        float yscale = -(float)historyRect[3]/(vars[i].max-vars[i].min);
        float yoffset = historyRect[1] + historyRect[3] - vars[i].min*yscale;
        stroke(darken(vars[i].thecolor));    strokeWeight(1.5);    noFill();
        beginShape();
        if (n0<n1) {
          for (int k=n0; k<=n1; k++) {
            vertex(historytime[k]*xscale+xoffset, constrain(varhistory[i][k]*yscale+yoffset,historyRect[1],historyRect[1]+historyRect[3]));
          }
        } else {
          for (int k=n0; k<length; k++) {
            vertex(historytime[k]*xscale+xoffset, constrain(varhistory[i][k]*yscale+yoffset,historyRect[1],historyRect[1]+historyRect[3]));
          }
          for (int k=0; k<=n1; k++) {
            vertex(historytime[k]*xscale+xoffset, constrain(varhistory[i][k]*yscale+yoffset,historyRect[1],historyRect[1]+historyRect[3]));
          }        
        }
        endShape();  
        textAlign(LEFT);    fill(vars[i].thecolor);    noStroke();
        textVAlign(" " + vars[i].shortname, historyRect[0]+historyRect[2], varhistory[i][n1]*yscale+yoffset, medfontsize, "middle");
      }
    } 
  }
  
  void draw_history() {
    if (showHistory) {
      int n1 = nhistend;
      int n0 = 0;
      if (wrapHistory) {
        n0 = nhistend+1;
        if (n0==length) {n0 = 0;} 
      }
      draw_history_setup(n0,n1);
      draw_history_content(n0,n1);  
    } 
  }
  
  
  void draw_other_things() {}
  
  void showDiagnostics() {}
  
  boolean over() {
    return (((mouseX >= networkRect[0]) && (mouseX <= networkRect[0]+networkRect[2])) && ((mouseY >= networkRect[1]) && (mouseY <= networkRect[1]+networkRect[3])));
  }

  
} // end of class dynamicalSystem




// -----------------------------------------------------------------------------------
//                  dynamical variables: state vars, fluxes, free params
// -----------------------------------------------------------------------------------



class flux { 
  int from, to; // indices of statevars; -1 = to/from nothingness
  float initial, current,scale,max,min;
  String name, shortname;
  color thecolor, currentcolor;
  float curvature, width, arrowheadpos;
  boolean display;
  
  flux(int N, int fillvalue) {
    display = false;
  }

  void define(String thename, String theshortname, int thefrom, int theto, color colr, float curv) {
    name = thename; shortname = theshortname;
    from = thefrom; to = theto;
    initial = 0;
    current = 0;
    thecolor = colr;
    curvature = curv/180*PI;
    arrowheadpos = 0.75;
    display = true;
    scale = 1;
    min = 0;
    max = big;
  }
  
}

class statevar {
  float initial,current,scale,max,min;
  String name, shortname;
  float xc,yc,radius;
  color thecolor, currentcolor;
  boolean constant,display;
  float sinkingRate = 0; // used by coupledBiophysicalModel; sinking is handled differently in the NNPZD and allo-upwelling1D models written before v1.3.re
  
  statevar(int N, int fillvalue) {
    constant = true;
    display = true;
  }
  
  void define(String thename, String theshortname, float initialval, color colr, float theX, float theY) {
    name = thename; shortname = theshortname;
    initial = initialval;
    current = initialval;
    constant = false;
    thecolor = colr;
    xc = theX; yc = theY;
    scale = 1;
    min = 0;
    max = big;
  }
    
}

class freeParam {
  float initial,min,max,current,previous;
  String name, shortname;
  color thecolor;
  slider controller;
  boolean controllerOn, frozen;
  
  freeParam(int N, int fillvalue) {
    frozen = false;
  }
  
  void define(String thename, String theshortname, float initialval, float minval, float maxval, color colr) {
    name = thename; shortname = theshortname;
    min = minval; max = maxval;
    controllerOn = false; // no slider by default
    initial = initialval;
    current = initialval;
    previous = initialval;
    thecolor = colr;
  }
  
  void addController(int[] screenRect) {
    controllerOn = true;
    controller = new slider(screenRect[0],screenRect[1],screenRect[2],screenRect[3],1);
    controller.addName(name);
    controller.highlightColor = thecolor;
  }
  
  void update() { // sync with controller & put result in _current_
    if (controllerOn) {
      if (frozen) {controller.freeze();} else {controller.unfreeze();}
      controller.setPos((current-min)/(max-min));
      controller.addDisplayVal(current,2);
      controller.update();
      current = controller.getPos() * (max-min) + min;
    }
  }
  
  void onMousePress() {}
  
}