//********************************************************************************* 
// Standard ORCHESTRA object class definitions
//
// Version 2021 (23 February)
//
// J.C.L. (Hans) Meeussen
// Nuclear Research and Consultancy Group, Petten, The Netherlands
// www.nrg.eu
// www.meeussen.nl/orchestra/objects2021.txt
// e-mail: meeussen@nrg.eu orchestra@meeussen.nl
// 
// This file contains the standard ORCHESTRA objects for implementing 
// chemical reaction and mass transport models.
// 
//
// Orchestra Objects (classes) are defined in terms of: 
// 1) Variables:                         
//	@Var: <name> <initial value>              
// 	@Var: pH 7
//
// 2) Mathematical expressions:        
//	@Calc:(1, "<variable> = <expression>") 
// 	@Calc:(1, "pH = -{H+.logact}")
//
// 3) Unknown - Equation pairs:              
//	@solve: unknown:(name:, pH, type:, lin) equation:(name:, H+.tot, tol: , 1e-15)
//
//
//
// The ENTITY, PHASE, and REACTION classes are the three basic building blocks
// for all chemical models.
//
//
// The definition of the phase and entity object classes can be found within the orchestra2.jar file.
//
// Changes in the 2021 version
// - The object classes @ini_phases() and @use_phases() that are used to initialise chemistry files were moved from the orchestra2021.jar file to this objectfile
// - This makes input files for  C++ and Java version identical.
// - (The C++ version only works with object files from this version 2021, or with the  @ini_phases() and @use_phases() provided in a separate file)
// - The expander (preprocesor) has a new keyword "@scan: <file>" that scans a file for class definitions without returning the text of the file (in contrast with @include: <file>) 
// - Filenames in @include: and @scan: statements can now contain object references e.g. @include: @chemistryfile() or @scan: @objectsfile()
// 
//  
// Changes in the 2020 version
// - Increased maximum ionic strenght to 20 moles/l inset_ion_activity_model class.
// - Adapted the species class to calculate activity for neutral species according to (logf = -0.1*I)
// - Added the SIT0 class to calculate SIT interaction for neutral species.
// - Added the @calculate_water_activity() class andadapted the species class 
//   this class can be called after including this file. 
//   Calculation of water activity is not switched on by default to remain compatible with older systems.
// - In this version the old (>20 year) "davies" and "species" objects were revised/simplified to allow easier replacement with
//   alternative ion activity models. (e.g. SIT).
//   In the new structure each ion(species) has its own activity coefficient <name>.logf, 
//   rather than the previously used common coefficient for all ions with a certain charge (f<charge>.log).
// - The SIT objects are currently present in a separate file (SITmodel2020.txt), that can be included in a chemistry file after this file.
// - The (obsolete) electroneutral diffusion objects have been removed, as these are superseded with the new
//   multicomponent / electrodiffusion classes.
//   
//
//
// Changes in the 2019 version
// - changed the settings of the equation in the nica_site object (solves problems with convergence)
// - introduced the <name>.<phase>.con variable
// - adsorption model objects are now defined in the same file as adsorption reactions
//
// Changes in the 2018 version
// - introduced the @solve: keyword as replacment for the @UnEq2: keyword for defining unknown-equation pairs
// - adsorption model specific code was removed from this file and added to adsmodels2018.txt
// - the nicamodel class was renamed to (more appropriate) donnanmodel class 
// - removed the maximum value constraint for surface site fractions (was 1)
// - the expander was extended with the possibility to evaluate expressions during the expanding process e.g: @evaluate:("2*3")
// - the expander does not require % characters in nested class definitions or sweeps anymore {%   {  }   %}
// 
//
//
// Changes in the 2014 version
// -  mineral object with variable tolerance that is adapted during iteration
// -  updated solid solution object
//   
// Changes in the 2013 version
//  - added a set of @logKreaction() objects which require 
//       a logK value as input rather than a k value. The logK reactions can be mixed with normal @reactions() in a single
//       input file to ensure backward compatibility.
//       The graphical chemistry editor understands both formats and will keep existing format by default.  
//  - The donnan model and nica sites were made electroneutral by adding dummy proton entities.
//    In this way the overall system remains electroneutral when surface area is changed. 
//    The oxide surfaces are are expressed in terms of neutral surface sites, but here the composition of 
//    the DDL is not explicitly taken into consideration. -->
//  - A mineral entity now has a .con variable to make it consistent with other entities, and to make it possible
//    to use it as a primary entity.
//  - A disperse() class was added, that calculates dispersion when combined with the convection object
//  - a balance_charge(cation,anion) object was added to use electroneutrality as constraint 
//
// changes in the 2011 version
//  - half life and decay objects added
//  - minerals are now phases as well as entities (made optional for efficiciency reasons)
//
//
// Changes in the 2008 version
//  - conversion to log activities
//  - changed all .act variables to .logact 
//  - removed all .sum variables (obsolete)
//  - changed activity correction factors to log equivalents (f-1.logact etc.)
//  - defined all "unkown" variables as global variables with default values
//  - changed solid solution iteration from log to linear
//  - added surface reaction with 6  components
//  - changed primary entity object for 7 parameters, (used for pH pe) so that it uses given default value for unknown
//  - updated and tested ddl objects (new version can handle charge unbalance in solution)
//**************************************************************************************

@class: ini_phases(){
 
  // original text from ini_phases.txt
 
  // we define the class phases, all phases are added to this class
   @Class: phases(){}

   // Each phase appends a bit of text to the phase class
   @Class: phase(phase){
      @Append: phases(){
         @Append: entity(name, phase){
            @Var:  <name>.<phase> 0
            @Var:  <name>.<phase>.con 0   //** Test
            @Calc:(1,"<name>.<phase>.con = 0")     //** Each entity gets a .con variable in each phase
            @Calc:(1,"<name>.<phase>     = 0")
         }
         @Append: link(name1, name2, coef1, coef2){
           @Calc:(3, "<name2>.<phase> = {<name2>.<phase>} + {<name1>.<phase>} * <coef2>")
         }
      }
   }

   @Class: phase(phase, parent, factor){
      @phase(<phase>)
      @link_phase(<phase>,<parent>,<factor>)
   }

   @Class: link_phase(phase, parent, factor){
      @Append: phases(){
         @Append: entity(name, phase){
            @Calc:(4,"<name>.<parent>     = {<name>.<parent>}     + {<name>.<phase>}     * <factor>")
            @Calc:(4,"<name>.<parent>.con = {<name>.<parent>.con} + {<name>.<phase>.con} * <factor>")  
         }
      }
   }
}

   // use_phases2.txt
@class: use_phases2(){
   @Class: entity(name, phase){
      @Var: <name>.logact 0
      @Var: <name>.con 0
      @Var: <name>.<phase> 0
      @Calc:(3,"<name>.<phase>.con = {<name>.con}")      //** test
   }

   @Class: entity(name, phase, factor){
      @entity(<name>, <phase>)
      @calc_conc(<name>, <factor>)
      @Calc:(2,"<name>.<phase>     = {<name>.con}")
   }

   @Class: entity(name, phase, ini, factor){
      @entity(<name>, <phase>, <factor>)
      @Calc:(1,"<name>.logact = <ini>")  
   }
   
   @Class: link(name1, name2, coef1, coef2){
      @Calc:(1, "<name1>.logact = {<name1>.logact} + {<name2>.logact} * <coef1>")
   }

   @Class: link(name1, name2, coef1){
      @link(<name1>, <name2>, <coef1>, <coef1>)
   }
   
   @phases()
}


@GlobalVar: T 298.15      // Temperature with default value

@Class: calc_conc(entity, factor){
   @Calc:(2,"<entity>.con = (10^{<entity>.logact}) * <factor>")
}

@Class: reaction(entity, k){
   @Var: <entity>.k   <k>
   @Calc:(1,"<entity>.logact = log10({<entity>.k})")
}

// The reaction classes with 1 - 10 reactants
@Class: reaction(entity, k, c1, n1){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
   @link(<entity>, <n8>, <c8>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
   @link(<entity>, <n8>, <c8>) 
   @link(<entity>, <n9>, <c9>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9, c10, n10){
   @reaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
   @link(<entity>, <n8>, <c8>) 
   @link(<entity>, <n9>, <c9>) 
   @link(<entity>, <n10>, <c10>) 
}


// The logk reaction allows to use logK values as input for reaction definitions 
// these can be mixed with ordinary reactions in input files

@Class: logKreaction(entity, logk){
   @Var: <entity>.logk   <logk>
   @Calc:(1,"<entity>.logact = {<entity>.logk}")
}

// The reaction classes with 1 - 8 reactants
@Class: logKreaction(entity, logk, c1, n1){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
}

@Class: logKreaction(entity, logk, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8){
   @logKreaction(<entity>, <logk>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
   @link(<entity>, <n8>, <c8>) 
}

@Class: logKreaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9){
   @logKreaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
   @link(<entity>, <n8>, <c8>) 
   @link(<entity>, <n9>, <c9>) 
}

@Class: logKreaction(entity, k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6, c7, n7, c8, n8, c9, n9, c10, n10){
   @logKreaction(<entity>, <k>)
   @link(<entity>, <n1>, <c1>) 
   @link(<entity>, <n2>, <c2>) 
   @link(<entity>, <n3>, <c3>) 
   @link(<entity>, <n4>, <c4>) 
   @link(<entity>, <n5>, <c5>) 
   @link(<entity>, <n6>, <c6>) 
   @link(<entity>, <n7>, <c7>) 
   @link(<entity>, <n8>, <c8>) 
   @link(<entity>, <n9>, <c9>) 
   @link(<entity>, <n10>, <c10>) 
}



//********************************************************************************* 
// A set of generic classes for constructing activity correction models
//
// Hans Meeussen  22 June 2020. 5 November 2020 
//********************************************************************************* 

@Class: set_ion_activity_model(name)
This object solves the unknown ionic strenght until estimated and calculatd ionic strength are equal.
It defines the ionic strength I as an unknown, and the difference between
the calculated and estimated ionic strength I.eq as the equation value which has to
become zero during the iteration.  
Hans Meeussen,  June 2002, revised March 2020
{
   @Var: I.eq 0                      // the equation that has to become zero
   @Calc: (3, "I.eq = I - {I.calc}") // the difference between calculated and estimated ionic strenght
   @<name>(0.1)
   @solve: unknown:(name:, I,  delta:, 1e-1, type:, lin, step: , .1, min:, 1e-9, max:, 20, iia:, true) equation:(name:, I.eq, tol:, 1e-4)
}


@Class: ion_activity_model(i)
This object calculates the activity correction factor I.act for ions charged 1.
The actual equation is provided in the class "@calculatef1()"
It also initializes the chargebalance and the calculated ionic strength (I.calc) variables, which are calculated in the species class.
It also calculates the electric conductivity EC.
The estimated/given ionic strength I is used as input.

Hans Meeussen,  June 2002, revised March 2020
{
   @globalvar:  I  <i>            // estimated ionic strength  
   @Var: I.logf1  0               // the log activity  coefficient for an ion with charge 1 
   @Var: I.calc   0               // I as calculated from actual speciation
   @Var: EC       0               // Electrical conductivity millimhos/cm
   @var: chargebalance 0          // chargebalance of dissolved species
   @Var: sumAqIons 0              // sum of ion concentrations used to estimate water activity

   @Calc: (1, "chargebalance = 0") // initialize chargebalance
   @Calc: (1, "I.calc        = 0") 
   @Calc: (1, "sumAqIons       = 0") 

   // Calculate the log activity  coefficient for an ion with charge 1 
   @Calc:(1, "I.logf1  = @calculatef1() "  )

   // Calculate electric conductivity from calculated I
   //EC (µS/cm)   =   6.2 10e4 × I (mol/L)   
   @Calc:(5, "EC      = 6.2e1 * I.calc")  



// we also define an allSpecies class to process all species 
@Class: allSpecies(){}


// This class redefines the 
// the species class to include the activity corrections

@Class: species(name, charge)
This new (March 2020) of the species class calculates its own activity coefficient, 
and contributions to charge balance and ionic strenght
{
   // we add this species to the list of species
   // so we can use this list to process all species
   @append: allSpecies(){
      @forEachSpecies(<name>, <charge>)
   }

   @entity(<name>, "@aqphase()")
   @Var: <name>.logf 0
   @Calc:(1,"<name>.logf   =  I.logf1 * (<charge>*<charge>) ")
//   @Calc:(1,"<name>.logf   =  if(<charge> == 0, (-0.1 * I),  I.logf1 * (<charge>*<charge>)) ")

   @Calc:(2,"<name>.con    = 10^({<name>.logact} + {<name>.logf}) ")
   @Calc:(2,"<name>.@aqphase() = {<name>.con}") 

   // add contribution of this species to total ionic strenght
   @Calc:(2,"I.calc  = I.calc + ({<name>.con} * (0.5*<charge>*<charge>))")

   // add contribution of this species to charge balance 
   @Calc:(2,"chargebalance = chargebalance + ({<name>.con} * (<charge>))")

   // add contribution of this species to sum of dissolved ions
   @Calc:(2,"sumAqIons = sumAqIons + {<name>.con}")
}


@Class: calculate_water_activity(){
   // calculation of water activity according to Phreeqc manual version 2 p17 https://pubs.usgs.gov/wri/1999/4259/report.pdf
   // water activity = 1 - 0.017 * sumAqIons
   @Var: H2O_activity_equation 0
   // we max the total mol/l of ions at 20 to prevent problems in predominace diagrams outside the water stability boundaries
   @Calc:(5,"H2O_activity_equation = H2O.logact - log10(1 - 0.017 * min(20, sumAqIons))" )
   @solve: unknown:(name:, H2O.logact,  delta:, 1e-3, type:, lin, step: , .1, reldelta:, 1e-3) equation:(name:, H2O_activity_equation, tol: , 1e-4 )
}

} 
//********************************************************************************* 




//********************************************************************************* 
// DAVIES model implementation
//
// Updated 22 June 2020, Hans Meeussen
//********************************************************************************* 
@Class: davies(){
   @set_ion_activity_model("davies")
}

@Class: davies(i){
   @Class: calculatef1(){ (0.51 * ( (sqrt(I)/(sqrt(I)+1))-(0.3*I))) }
   @ion_activity_model(<i>)
}
//********************************************************************************* 



//********************************************************************************* 
// Complete ORCHESTRA SIT model implementation, plus set of SIT parameters
//
// To use the SIT model, include this file after the standard objects file.
// an select SIT activity correction model in the GUI.
//
// Hans Meeussen,   6 October 2020, tested against PHREEQC, 23 Februari 2021 
//********************************************************************************* 

@Class: sitmodel(){
   @set_ion_activity_model("sitmodel")
}

@Class: sitmodel(i){

@Class: calculatef1(){ (0.51*sqrt(I)) / (1+1.5*sqrt(I))}

@ion_activity_model(<i>)

@class: SIT(ion1, ion2, coef){
   @var: <ion1>.logf     0
   @var: <ion2>.logf     0
   @var: <ion1>.logact -30
   @var: <ion2>.logact -30
   @var: <ion1>.logfsum  0
   @var: <ion2>.logfsum  0

   // we calculate the logsum from all the interactions and concentrations
   // only contributions for ions with a log activity >-4

   @Calc: (1,"{<ion1>.logfsum} = {<ion1>.logfsum} + if({<ion2>.logact}>-4,(<coef> * 10^({<ion2>.logact} + {<ion2>.logf})), 0)")
   @Calc: (1,"{<ion2>.logfsum} = {<ion2>.logfsum} + if({<ion1>.logact}>-4,(<coef> * 10^({<ion1>.logact} + {<ion1>.logf})), 0)")
}

// called for each species
@class: SIT_update(name, charge){
   @var: <name>.logfsum 0

//   if charge == 0, we could use a specific equation here, but this is currently not used
//   @Calc:(1, "<name>.logf = if(<charge> == 0 , (-.01*I) - {<name>.logfsum}, (<charge>^2) * @calculatef1() - {<name>.logfsum})")

// apparently PHREEQC uses zero values for logf of neutral species, so we use the equation below
//  @Calc:(1, "<name>.logf =  (<charge>^2) * @calculatef1() - {<name>.logfsum}")
   @Calc:(1, "<name>.logf =  (<charge>^2) * I.logf1 - {<name>.logfsum}")
   
   @Calc:(1, "<name>.logfsum = 0")
}

// an SIT interaction coefficient for neutral species with background electrolyte ions. 
@class: SIT0(neutral_species, cation, anion, coef){
   @var: <neutral_species>.logf     0
   @var: <neutral_species>.logfsum  0

   // We use the minimum concentration of the given anion/cation combination, so we do not wrongly add up in case of mixed background electrolyte (e.g. N+, K+, Cl-)
   // if this concentration is <1e-4 we ignore it
   @Calc: (1,"{<neutral_species>.logfsum} = {<neutral_species>.logfsum} + if(min({<cation>.logact}, {<anion>.logact}) >-4, <coef> * 10^min(({<cation>.logact} + {<cation>.logf}), ({<anion>.logact}  + {<anion>.logf})) ,0)")
}


// this class needs to be called from the chemistry file after the species definition
// the ORCHESTRA GUI automatically adds this to chemistry files when the SIT model is selected
@Class: sitparameters(){
 
   // we define the forEachSpecies task that is called for each species
   // when calling the allSpecies object below
   @class:  forEachSpecies(name, charge){
      @SIT_update(<name>, <charge>)
   }

   // calculate the SIT sum factors using current activity coefficients
   @sitparameters2()

   // update logf factors with sum factors 
   @allSpecies()

   
   // repeat this a number of times, this converges very quickly 
   @sitparameters2()
   @allSpecies()

   @sitparameters2()
   @allSpecies()

   @sitparameters2()
   @allSpecies()

}

@Class: sitparameters2(){
  // this is a default empty set of sitparameters
  // a full set can be found and used by scanning the SIT database file after inlcluding this objects file
  //@scan: sitparameters.txt
  // we define the class sitparameter here, so this is automatically done when we 
  // use the sitmodel class

@SIT(CO3-2,Na+,-0.08)
@SIT(Ca+2,Cl-,0.14)
@SIT(Ca+2,NO3-,0.02)
@SIT(Cl-,Na+,0.03)
@SIT(H+,Cl-,0.12)
@SIT(H+,NO3-,0.07)
}

}

// The default aqueous phase is "diss", however this can be redefined if necessary. 
@Class: aqphase(){diss}



@Class: mineral(name).
 
This object represents an automatically precipitating/dissolving mineral.
A mineral is an entity in the "min" phase, which can be combined with a standard (formation)reaction object. 
The ".logact" variable of a mineral entity is its log IAP + log K value.

Because, in contrast with ordinary entities, the amount of this mineral entity cannot be calculated directly from its activity we 
introduce an extra unknown and equation that is solved by the standard iteration procedure.
The amount of mineral can now be calculated from the estimated unknown.
If the unknown is negative, the mineral is not present and the unknown represents its log undersaturation.
If the unknown is positive it directly represens the amount of mineral present.

The equation that has to become zero to reach convergence:
unknown <  0  equation = log (IAP*K) - unknown   // difference between unknown and log undersaturation  (should go to zero)
unknown >= 0  equation = log (IAP*K)             // log saturation index                                (should go to zero)

In this way we can implement a precipitating mineral object calss with just 10 statements 

May 2014
Hans Meeussen.

{
     //   @phase(<name>, min, 1)
	@entity(<name>, min) 
	@Var: <name>.un   -1e-3
	@Var: <name>.eq    0
	@Var: <name>.si    0
	@Var: minTol      1e-3  // minTol goes to zero during the iteration

	@Calc:(2,"<name>.si  = {<name>.logact}") 
	@Calc:(2,"<name>.con = if({<name>.un} >= 0,  ({<name>.un}),   0)  ")
	@Calc:(2,"<name>.min = {<name>.con}")

	@Calc:(2,"<name>.eq  = if({<name>.un} < 0,     ({<name>.si} - {<name>.un}), {<name>.si}- {<name>.un} * minTol )")
	@Solve: unknown:(name:, <name>.un,  delta:, 1e-6, type:, lin, step: , .1, reldelta:, 1e-3) equation:(name:, <name>.eq, tol: , 1e-8, si:, <name>.si )
}




@Class: solid_solution(name, parent_phase)
Solid solutions in are solved in ORCHESTRA by introducing an extra unknown and equation.
When negative the unknown represents the undersaturation, and indicates that the amount of solid solutio is zero.
The accompanying equation then is the difference between unknown and (log) undersaturation.
When positive the unknown represents the amount of solid solution. 
The accomponying equation then is the difference between 1 and the sum of fractions

{
	@GlobalVar:  <name>.un        -1  // the unknown that is used in the iteration
	@GlobalVar:  <name>.eq         0  // the calculated equation that should become zero during iteration
	@Var:  <name>.est_sum          0  // the estimated sum (is calculated from unknown)

	// a solid sulution is a phase in which entities can be present
	@link_phase(<name>, <parent_phase>, <name>.est_sum)

	// a solid solution is also an entity (in its own phase)
	//(so has an activity and mass balance and can be used in reactions)  
     // logact = 0 
     @Globalvar: <name>.logact 0
     @entity(<name>, <name>, 0, 0)

	// calculate estimated amount of solid solution (est_sum) directly from the unknown
	@Calc:(1,"<name>.est_sum = if({<name>.un} >= 0,   {<name>.un},   0)  ")

	// calculate equation (error that should become zero) from unknown and calculated sum of fractions
     // this implies that the sum of fractions needs to become 1 when solid solution exists.
     // if solid solution does not exist  (unknown<0), the equation is undersaturion.

	// if solid solution exists,                  (unknown>0)  estimated sum should be equal to calculated sum 
	// if solid solution is completely dissolved, (unknown<=0) calculated sum should be zero
     // the variable <name>.<name> contains the sum of fractions that should go to 1.

      @Calc:(5, "<name>.eq  = if(0 >  {<name>.un},   ( log({<name>.<name>})/1000 - {<name>.un}), 1-{<name>.<name>})")

	// iterate for a value of the unknown that results in an equation value of zero.
      @solve: unknown:(name:, <name>.un, delta:, 1e-3, type:, lin, step: , .1 ) equation:(name:, <name>.eq, tol: , 1e-6)

}



@Class: e_layer (name, fase)
An entity that represents an electrostatic layer with 
<name>.logact         = log boltzmann factor 
<name>.<phase>        = charge balance. 
By placing this entity in an e fase, its sum is directly converted to the units Coulomb/m2

{
	@entity (<name>, <fase>, 0)
	@Var: <name>.psi .1
   	@GlobalVar: T 298.15      // Temperature with default value

	//* calculate psi from boltzman factor name.act
	@Calc:(1,"<name>.psi = (-{<name>.logact}* log(10))/(96484.56/(8.31441*T))") 
}


@Class: ddl_notyetcompatible(name, e_layer)   
{ 
    // Object to calculate diffuse double layer charge from given potential of electrostatic layer
    @Var: <name>.ddl 0.0  
    @Var: <name>.PF_RT 1
    @GlobalVar: T 298.15          // Temperature with default value
    @Var: R    8.31441
    @Var: F    96484.56
    @Var: eps0 8.85419e-12        // dielectric constant vacuum Sposito p 229
    @Var: eps1 78                 // dielectric constant water Hiemstra p 46

    @Calc:(1,"<name>.PF_RT = {<e_layer>.psi} * (F/(R*T)) ")  

    @Var: <name>.sum_of_charge  0

    @Calc:(4,"<name>.ddl        = if({<e_layer>.psi} > 0, {<name>.ddl}, -{<name>.ddl})")
    @Calc:(4,"<name>.ddl        = if({<name>.sum_of_charge} > 0, -sqrt(2000*eps0*eps1*R*T) * sqrt({<name>.sum_of_charge}), 0 ))")

      
    // if the charge balance of the solution is not zero (which can happen during iterations, or partial system definition)  
    // we assume that the charge unbalance is due to a missing -1 or +1 species and add its contribution to the sum 
    @Calc: (4, "<name>.sum_of_charge  = <name>.sum_of_charge  +  if (chargebalance>0,  chargebalance * (exp(<name>.PF_RT) -1), -chargebalance * (exp(-<name>.PF_RT) -1) )"

    // we use the total concentrations of dissolved species for each distinct charge, as calculated by the ionic strength object
    @Calc:(4,"<name>.sum_of_charge  = {f1.sum} * (exp(-1 * <name>.PF_RT) -1) + {f-1.sum} * (exp(<name>.PF_RT) -1) + {f2.sum} * (exp(-2 * <name>.PF_RT) -1) + {f-2.sum} * (exp(2 * <name>.PF_RT) -1) + {f3.sum} * (exp(-3 * <name>.PF_RT) -1) + {f-3.sum} * (exp(3 * <name>.PF_RT) -1)")
}


@Class: ddl (name, e_layer){
    // Object to calculate diffuse double layer charge from given potential of electrostatic layer
    // This version of ddl object was used in first Orchestra versions and assumes that the ionic strenght is dominated by symmetric electrolyte monovalent ions
    @Var: T    298.15          
    @Var: R    8.31441
    @Var: F    96484.56
     
    @Var: <name>.ddl 0  
    @Calc:(1,"<name>.ddl  =  -0.0587 * sqrt( I * (exp(-1 * {<e_layer>.psi} * (F/(R*T)) ) -1) + I * ( 1 * exp({<e_layer>.psi} * (F/(R*T)) ) -1))   )")    
    @Calc:(1,"<name>.ddl  =  if({<e_layer>.psi} > 0, {<name>.ddl}, -{<name>.ddl})")
 }


@Class: ddl_dzombak (name, e_layer){
        @Var: F    96484.56
      // This version of ddl model is used by the Generalized two layer model of Dzombak and Morel.
	@Var: <name>.ddl 0.0  
	@Calc:(4,"<name>.ddl = -2.5*sqrt(I)*{<e_layer>.psi} ")  //dzombak & morel p12
}



//**********************************************************************************
// The transport object classes
// *********************************************************************************

//------------------------------------------------------------------------------------
// Transport objects translated to new expressions (9 September 2002)

@Class: update_mass (name) 
Updates the mass of a component in dynamic systems (transport or kinetics). 
It is normally used within the file update_mass.inp.
It assumes that the phases "tot" and "d" are used to represent total and 
delta (per time unit) amounts of component masses.
<br> Hans Meeussen
{  @globalvar: <name>.tot	0
   @globalvar: <name>.d	0
   @Var: dt       		1

   // Add the mass change <name>.d (mol/s) multiplied by the time step dt(s) 
   // to the total amount in the cell and reset name.d.sum to zero.
   @Calc:(1, "<name>.tot = {<name>.tot} + {<name>.d} * dt")
   @Calc:(1, "<name>.d = 0")
}

@Class: convection ()
This object contains the general part of convection and defines the 
amount of water transported between cell 1 and cell 2.
It should be used in convection.inp before the component specific objects.
{  @Var: dwater  0    
   @Var: 1.J     0
    
   //The amount of water moving from cell 1 to cell 2 equals the 
   //flow rate (l/s). 
   @Calc:(1, "dwater=1.J")
}

@Class: convec (name) 
This object contains the component specific part of convection and defines the 
amount of component (mol/s) transported between cell 1 and cell 2.
{  @globalvar: 1.<name>.diss 0
   @globalvar: 2.<name>.diss 0
   @Var: 1.<name>.d 0
   @Var: 2.<name>.d 0
   @Var: dmass 0 

   //Calculate the transported mass between cell 1 and cell 2 
   //depending on the direction of the flow
   @Calc:(1, "dmass = if (dwater>0 , {1.<name>.diss} * dwater, {2.<name>.diss} * dwater)")
   @Calc:(1,"1.<name>.d = {1.<name>.d} - dmass")
   @Calc:(1,"2.<name>.d = {2.<name>.d} + dmass")
}

@Class: diffuse (name, D) {
   @Class: diffuse (<name>, <D>, diss) 
}

@Class: diffuse (name, D, phase) 
This object is used to calculate diffusion of a solute driven by its 
concentration gradient according to Fick's law. It uses the
concentrations in both cells, the distance between the cells (dx)
in meter, and the surface area between the cells (A) in m2.
<br>Wendy van Beinum, Hans Meeussen
{  @Var: 1.<name>.<phase>	0
   @Var: 2.<name>.<phase>	0
   @Var: 1.<name>.d     0
   @Var: 2.<name>.d     0
   @Var: <name>.D  <D>


   //Calculate concentration gradient in (mol/m3)/m
   //dc=1000(C2-C1)/dx
   @Var: <name>.dc 0
   @Var: 1.dx      1
   @Calc:(1, "<name>.dc = 1000*({2.<name>.<phase>} - {1.<name>.<phase>})/1.dx")

   //Calculate the transported amount in mol/s
   //J= -D*dC*A 
   @Var: <name>.J  0
   @Var: 2.A       0
   @Calc:(1, "<name>.J=(-{<name>.D} * {<name>.dc})*2.A")

   //Add and substract the transported amount from the mass changes in both cells
   @Calc:(1, "1.<name>.d = {1.<name>.d} - {<name>.J}")
   @Calc:(1, "2.<name>.d = {2.<name>.d} + {<name>.J}")
}


@Class: disperse (name, dispersivity, phase) 
This object is used to calculate dispersion of a solute driven by its 
concentration gradient, (material) dispersivity and pore-water velocity.
It uses the concentrations in both cells, the distance between the cells (dx)in meter, 
the surface area between the cells (A) in m2, and the pore-water velocity (A*(J/porosity)) as input.
J = in liters/s
J/A = liter per s per m2
A is total contact area
D = m
porosity = water filled porosity


//<br> Gijsbert Cirkel, Hans Meeussen   25-09-2012
{  @Var: 1.<name>.<phase>	0
   @Var: 2.<name>.<phase>	0
   @Var: 1.<name>.d     0
   @Var: 2.<name>.d     0
   @Var: <name>.D  <dispersivity>
   @Var: porosity   0.38   
   @Var: 1.porosity   1
   @Var: 2.porosity   1

//**** 
   The dx and porosity calculations are generic for all substances, 
   and could bedefined outside this substance specific object.
   However, defining them here will not affect performance 
   as the ORCHESTRA optimizer will remove the redundant calculations.
*******//

   @Var: 1.dx         0.001   //distance between cell 1 and 2 (m) 
   @Var: dx 1 // the local dx variable with default value
   @Calc:(1, "dx = 1.dx") // copy dx from cell

   @Calc:(1, "porosity = (1.porosity+2.porosity)/2")

//****

   //Calculate concentration gradient in (mol/m3)/m
   //dc=1000(C2-C1)/dx
   @Var: <name>.dc 0
   @Calc:(1, "<name>.dc = 1000*({2.<name>.<phase>} - {1.<name>.<phase>})")

   //Calculate the transported amount in mol/s
   //J= -D*dC*A*pore-water velocity 
   @Var: 1.A       1
   @Var: <name>.J       0

   // <name>.J = delta mass for this component, 
   // dwater =  amount of water in liter/sec
   // dwater/(A*porosity)  = pore-water velocity 
   @Calc:(1, "<name>.J=(-porosity * {<name>.D} * {<name>.dc}/dx) * 1.A * (abs(dwater)/(1.A * porosity))")  //porosity added GC 20-09-2012
//   @Calc:(1, "<name>.J=(-{<name>.D} * {<name>.dc}/dx) * abs(dwater)")  //simplified version

   //Add and substract the transported amount from the mass changes per time in both cells
   @Calc:(1, "1.<name>.d = {1.<name>.d} - {<name>.J}")
   @Calc:(1, "2.<name>.d = {2.<name>.d} + {<name>.J}")
}

// The user interface objects (used by the graphical chemistry editor)

@Class: primary_entity(name, unknownvalue, phase, equationvalue){
     @Var: tolerance 1e-13
     @globalvar: <name>.logact <unknownvalue>  

     @globalvar: <name>.<phase> <equationvalue>     // define this variable as global (node) variables + default value
     @solve: unknown:(name:, <name>.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, -15) equation:(name:, <name>.<phase>, tol: , tolerance, default:, <equationvalue>) )
}

@Class: primary_entity(name, unknownvalue){
     @globalvar: <name>.logact <unknownvalue>  
}

@Class: primary_entity(name, unknown, unknownvalue, type, step, phase, equationvalue){
     @Var: tolerance 1e-13
     @GlobalVar: <name>.<phase> <equationvalue>     // define this variable as global (node) variables + default value
     @solve: unknown:(name:, <unknown>, delta:,  1e-6, type:, <type>, step:, <step>, default:, <unknownvalue>) equation:(name:, <name>.<phase>, default:, <equationvalue>, tol:, tolerance)
}

@Class: primary_entity(name, unknown, unknownvalue){
     @globalvar: <unknown> <unknownvalue>  
}


@class: primary_entity_synonym:(synonym, name){
   @GlobalVar: <name>.tot 0
   @GlobalVar: <name>.diss 0
   @GlobalVar: <name>.ads 0
   @GlobalVar: <name>.min 0
   @GlobalVar: <name>.gas 0
   @Synonym: <synonym>.tot  <name>.tot
   @Synonym: <synonym>.diss <name>.diss
   @Synonym: <synonym>.ads  <name>.ads
   @Synonym: <synonym>.min  <name>.min
   @Synonym: <synonym>.gas  <name>.gas
}


// An imineral is just a mineral that is selected as a primary entity
// this creates a fixed activity 
@class: imineral(name){
    @entity(<name>, min, 1)
}

@Class: gas(name){
   @entity(<name>, gas, 1)
}



@Class: donnanmodel(name, parentphase, concentration, donnanVolume){
   // define the donnan surface phase
   @phase(<name>, <parentphase>, <concentration>)
   
   //--------------------------------------------------------------------------------------
   // Add a donnan phase to this surface, which consist of a phase, linked to the 
   // surface via volume/kg. The donnan phase is also an entity (logact = boltzman factor, sum: <name>_don.<name> = charge balance) 
   // unknown is activity, equation is charge balance at surface = 0 
   //--------------------------------------------------------------------------------------
   @Var: <name>donvol 1.5  // The Donnan volume + default value
   @Calc:(1,"<name>donvol = <donnanVolume>")
   @phase(<name>_don, <name>, <name>donvol)

   @Var: <name>_don.logact -1 
   @entity(<name>_don, <name>_don, 0)    
   @solve: unknown:(name:, <name>_don.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, 0, iia:, true) equation:(name:, <name>_don.<name>, tol: , 1e-8)
   //--------------------------------------------------------------------------------------------
}



@Class: donnanmodel(name, parentphase, concentration, donnanVolume, charge){
   @donnanmodel(<name>, <parentphase>, <concentration>, <donnanVolume>)
   @calc:(3,"<name>_don.<name> = {<name>_don.<name>} + <charge>")


   // by adding a dummy phase and dummy proton entity to this surface e make that the net surface charge zero.
   // otherwise surface would be negative and consumes protons, or cations upon addition to solution  
      
   //    @Calc:(1,"<name>donvol = <donnanVolume>")
   //    @phase(<name>_dummy, <name>, <name>donvol)

   //    @entity(<name>dummyproton, <name>_dummy, 1)
   //    @link(<name>dummyproton, H+, 0, <charge> ) 
}


@Class: nicamodel(name, parentphase, concentration, donnanVolume){
   @donnanmodel(<name>, <parentphase>, <concentration>, <donnanVolume>)
}

@Class: nicamodel(name, parentphase, concentration, donnanVolume, charge){
   @donnanmodel(<name>, <parentphase>, <concentration>, <donnanVolume>, <charge>)
}


//----------------------------------------------------------------------------------------------
// A NICA site can be represented by a standard phase plus a standard entity.
// The phase is linked to the total particle surface phase, by the site density sd (moles/kg )
//----------------------------------------------------------------------------------------------
@Class: nicasite(name, surface, donnan, p, sd)
{
	// A nica site is a phase that is linked to its parent surface via its site densitiy
	@phase(<name>, <surface>, <sd>)

	// it is also an entity in its own phase
	@entity(<name>, <name>, 0)

        // by adding a dummy proton entity to this site we make that the net surface charge of nica site = zero.
        // otherwise adding nica-donnan surface is negative and consumes protons, or cations upon addition to solution  
	@entity(<name>dummyproton, <name>, 1)
	@link(<name>dummyproton, H+, 0, -1 ) 

	// The amount of charge that is carried by the empty sites is added to the donnan mass (charge) balance 
	@Calc:(3,"<donnan>.<surface> = {<donnan>.<surface>} - <sd>")

	// The value of the unknown C is used to calculate C^(p-1)/(1+C^p) which is the .act of the NICA site.
	// The following lines are the following formula:  act = (unknown^(p-1))/(1+unknown^p)
	
	@Var: <name>.unknown 1
	@Var: <name>.equation 0

	@Calc:(1,"<name>.logact = log10(({<name>.unknown}^(<p>-1))/(1+{<name>.unknown}^<p>))")

	// Subtract the calculated sum(K^n* C^n) from the estimated sum, result should become zero
	// equation = calculated sum ((nic.act*C)/nic.act = C) - estimated sum (unknown)

	// equation is calculated logC - estimated logC, should become zero
	@Calc:(3,"<name>.equation = log10({<name>.<name>} / (10^{<name>.logact})) - log10({<name>.unknown})")

	@solve: unknown:(name:, <name>.unknown,  delta:, 1e-3, type:, log, iia:, true) equation:(name:, <name>.equation, tol: , 1e-6)
}



@Class: nicaspecies(name, site, ion, n, nH, logK)
{
   @entity(<name>, <site>, "<n>/<nH>")
   @Calc:(1, "<name>.logact = <logK>*<n>")

   @link(<name>,    <site>,   1, "(<nH>/<n>)")
   @link(<name>,    <ion>,  <n>,    1 )
}


@Class: donnanspecies(name, phase, ion, charge)
{	
   @entity(<name>, <phase>, 1)

   // here we use the ion concentration as input rather than the ion activity
   // so we multiply activity times activity coefficient
   @calc:(1,"<name>.logact = {<ion>.logf})") 
   @link(<name>, <ion>, 1)
   @link(<name>, <phase>, <charge>)
}

@class: adsmodel(name, parent_phase ,concentration, type){
   @<type>(<name>, <parent_phase> ,<concentration>)
}

@class: surfsite(model, name, density, coef){
   @phase(<model>_<name>, <model>, <density>)
   @entity(<model>_<name>, <model>_<name>, <coef>) 
   @solve: unknown:(name:, <model>_<name>.logact,  delta:, 1e-6, type:, lin, step: , .1, default:, 1) equation:(name:, <model>_<name>.<model>_<name>, tol: , 1e-4, default:, 1)
}

@class: surfsite(model, name, density, coef, charge, plane){
   @phase(<model>_<name>, <model>, <density>)
   @entity(<model>_<name>, <model>_<name>, <coef>) 
   @link(<model>_<name>, <model>_<plane>, 0, <charge>)
   @solve: unknown:(name:, <model>_<name>.logact,  delta:, 1e-6, type:, lin, step: , .1, default:, -1) equation:(name:, <model>_<name>.<model>_<name>,  tol: , 1e-4, default:, 1)
}

@class: surfsite(model, name, density, coef, charge, plane, charge2, plane2){
   @phase(<model>_<name>, <model>, <density>)
   @entity(<model>_<name>, <model>_<name>, <coef>) 
   @link(<model>_<name>, <model>_<plane>, 0, <charge>)
   @link(<model>_<name>, <model>_<plane2>, 0, <charge2>)
   @solve: unknown:(name:, <model>_<name>.logact,  delta:, 1e-6, type:, lin, step: , .1, default:, -1) equation:(name:, <model>_<name>.<model>_<name>,  tol: , 1e-4, default:, 1)
}


@class: surfspecies(model, site, name, coef){
   @entity(<model>_<name>, <model>_<site>, <coef>)
}

@class: surfreaction(name,  k, c1, n1){
   @reaction(<name>,  <k>, <c1>, <n1>)
}

@class: surfreaction(name,  k, c1, n1, c2, n2, ){
   @reaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>)
}

@class: surfreaction(name,  k, c1, n1, c2, n2, c3, n3){
   @reaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>)
}

@class: surfreaction(name,  k, c1, n1, c2, n2, c3, n3, c4, n4 ){
   @reaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>, <c4>, <n4>)
}

@class: surfreaction(name,  k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5 ){
   @reaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>, <c4>, <n4>, <c5>, <n5>)
}

@class: surfreaction(name,  k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6 ){
   @reaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>, <c4>, <n4>, <c5>, <n5>, <c6>, <n6>)
}


@class: logKsurfreaction(name,  k, c1, n1){
   @logKreaction(<name>,  <k>, <c1>, <n1>)
}

@class: logKsurfreaction(name,  k, c1, n1, c2, n2, ){
   @logKreaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>)
}

@class: logKsurfreaction(name,  k, c1, n1, c2, n2, c3, n3){
   @logKreaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>)
}

@class: logKsurfreaction(name,  k, c1, n1, c2, n2, c3, n3, c4, n4 ){
   @logKreaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>, <c4>, <n4>)
}

@class: logKsurfreaction(name,  k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5 ){
   @logKreaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>, <c4>, <n4>, <c5>, <n5>)
}

@class: logKsurfreaction(name,  k, c1, n1, c2, n2, c3, n3, c4, n4, c5, n5, c6, n6 ){
   @logKreaction(<name>,  <k>, <c1>, <n1>, <c2>, <n2>, <c3>, <n3>, <c4>, <n4>, <c5>, <n5>, <c6>, <n6>)
}





//**
// This object creates a kd entity and phase for a component that contains 
// an amount of mass equal to total dissolved concentration * kd
//
// usage in chemical input file: 
// @kd(K+, 4.8)
//
@class: kd(name, kd){
   @phase(<name>_kd, ads, "(<kd>*<name>.diss)")   // concentration in ads phase is dissolved concentration * kd
   @entity(<name>_kd,<name>_kd, 1)
   @link(<name>_kd, <name>, 0, 1)            
}


// This object calculates the amount of adsorbed protons necessary to agree with measured ANC
// usage in chemical input file: 
// @ANC("4e-3*pH")
//
@class: ANC(formula){
   // here we need a formula that calculates the total amount of adsorbed protons from the pH.
   // subtract the actual adsorbed amount of protons from this, so the effective total equals ANC
   @phase(ANC, ads, <formula> - H+.ads.sum")
   @entity(ANC, ANC, 1)
   @link(<name>_kd, <name>, 0, 1)            
}


@Class: balance_charge(cation,anion)
// usage: @balance_charge(Na+, Cl-)
// This object balances charge by adding either the amount of a cation or an anion, whichever is appropriate.
// The charge_min_act variable represents the lower limit value for the ion that is not corrected. 
// Make sure that for each of these ions a constant activity is chosen in the primary entities settings
// so no iteration on their individual mass balances is carried out.
// Hans Meeussen. July 2012
{
   @Var: chargebalance 0  
   @Var: charge_unknown -1e-3 // the unknown that is varied until chargebalance is zero.
   @Var: charge_min_act 1e-3  // the concentration of the cat/anion that is not corrected.

   @Calc: (1, "<cation>.logact = if(charge_unknown>0, log10(charge_min_act + charge_unknown), log10(charge_min_act))")
   @Calc: (1, "<anion>.logact  = if(charge_unknown<0, log10(charge_min_act - charge_unknown), log10(charge_min_act))")
   @solve: unknown:(name:, charge_unknown,  delta:, 1e-3, type:, lin, step: , .1) equation:(name:, chargebalance, tol: , 1e-6)
}


@Class: balance_charge(unknown) 
{
   @GlobalVar: chargebalance 0  
   @Calc: (1, "chargebalance = 0")   
   @solve: unknown:(name:, <unknown>,  delta:, 1e-6, type:, lin, step: , .1,iia:, true) equation:(name:, chargebalance, tol: , 1e-10, default:, 0)
}




// the radioactive decay objects
// Hans Meeussen, NRG, September, 2010

@class: radionuclide(name, value){
   @var: <name>.halflife   <value>  // half life in seconds
   @globalvar: <name>.tot  0        // The total amount of mother isotope
   @globalvar: <name>.d    0        // The delta amount of mother isotope
   @calc:(1,   "<name>.d     = {<name>.d}     - (log(2)/{<name>.halflife} * if({<name>.tot} > 1e-20, {<name>.tot}, 0))")
}

@class: halflife(name, value){
   @radionuclide(<name>, <value>)
}

@class: decay(mother, frac, daughter)
// mother component, fraction daughter component
{
   @globalvar: <daughter>.d 0   // The delta amount of daughter product

   // here we limit the the decay at a minimum total concentration of 1e-20 M to prevent zero concentrations
   @calc:(1,   "<daughter>.d   = {<daughter>.d}   + <frac> * (log(2)/{<mother>.halflife} * if({<mother>.tot} > 1e-20, {<mother>.tot}, 0))")
}

// Merge isotope concentrations into total element concentrations
@class: prechem_perphase(phase, isotope, component){
  @globalVar: <isotope>.<phase>  0
  @Calc: (1, "<component>.<phase> = {<component>.<phase>} + {<isotope>.<phase>}") 
}

// Split element concentrations back into individual isotope concentrations
@class: postchem_perphase(phase, isotope, component){
  @globalVar:   <isotope>.<phase>  0
  @globalVar: <component>.<phase>  0
  @Calc: (1, "<isotope>.<phase> = {<component>.<phase>}  * {<isotope>.frac}") 
}