//********************************************************************************* 
// Standard ORCHESTRA object definitions
//
// Version 2008
//
// Hans Meeussen
// Energy Research Centre of the Netherlands
// www.ecn.nl
// meeussen@ecn.nl
// 
// This file contains a standard set of ORCHESTRA object definitions that represent 
// (building blocks for) chemical reaction and mass transport models. 
// This set is not complete, but illustrates how the framework of basic 
// ORCHESTRA objects; ENTITY, PHASE, and REACTION, can be used to implement
// chemical models in a compact, consistent way.
//
//
// Orchestra Objects (classes) can be defined in terms of: 
// 1) Variables:                         
//	@Var:    <name> <initial value>              
// 	@Var: pH 7
//
// 2) Mathematical expressions:        
//	@Calc:(1, "<variable> = <expression>") 
// 	@Calc:(1, "pH = 10^-H+.act")
//
// 3) Unknown - equation pairs:              
//	@Uneq2: unknown:(name:, pH, type:, lin) equation:(name:, H+.tot, tol: , 1e-15)
//
// The definition of the phase and entity object classes can be found in the orchestra2.jar file
//
// Changes in the 2008 version
//  - conversion to log activities
//  - changed all .act variables to .logact 
//  - removed all .sum variables (not used anymore)
//  - 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)
//**************************************************************************************

@logactivities: // indicates that calculator will use log activities

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

@Class: calc_conc_with_davies(name, charge, phase){
	@Calc:(2,"<name>.con =	(10^({<name>.logact} + {f<charge>.log}))")
	@Calc:(2,"f<charge>.sum = {f<charge>.sum} + {<name>.con}")
	@Calc:(2,"<name>.<phase> = {<name>.con}") 
}

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

// The reaction classes with 1 - 5 reactants
@Class: reaction(entity, k, c1, n1){
	@Calc:(1,"<entity>.si2 = abs(<c1>)")
	@reaction(<entity>, <k>)
	@link(<entity>, <n1>, <c1>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>)")
	@reaction(<entity>, <k>)
	@link(<entity>, <n1>, <c1>) 
	@link(<entity>, <n2>, <c2>) 
}

@Class: reaction(entity, k, c1, n1, c2, n2, c3, n3){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>) + abs(<c3>)")
	@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){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>) + abs(<c3>) + abs(<c4>)")
	@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){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>) + abs(<c3>) + abs(<c4>) + abs(<c5>)")
	@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){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>) + abs(<c3>) + abs(<c4>) + abs(<c5>) + abs(<c6>)")
	@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){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>) + abs(<c3>) + abs(<c4>) + abs(<c5>) + abs(<c6>) + abs(<c7>)")
	@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){
	@Calc:(1,"<entity>.si2 = abs(<c1>) + abs(<c2>) + abs(<c3>) + abs(<c4>) + abs(<c5>) + abs(<c6>) + abs(<c7>) + abs(<c8>)")
	@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: davies()
This object creates a chemical system with a calculated ionic strength.
It defines the ionic strength I as an unknown, and the difference between
the calculated and estimated ionic strength I.sum as the equation value which has to
become zero during the iteration.  
{
	@davies(0.1)
	@Uneq2: unknown:(name:, I,  delta:, 1e-6, type:, lin, step: , .1, min:, 0, max:, 10, iia:, true) equation:(name:, I.sum, tol: , 1e-6)
}


@Class: f(charge) 
This object represents the ion activity correction factor for
ions, according to the Davies model. 
The "activity" (f.act) of this object represents the activity factor,
the "sum" f<charge>.sum the total amount of ions with this charge. 
<br>Hans Meeussen
{
   //***** act = activity coefficient, sum = sum * Z^2 of ions with charge z
	@Var: f<charge>.log 	0
	@Var: f<charge>.sum	0
	@Calc:(1,"f<charge>.sum = 0")
	@Calc:(1,"f<charge>.log = log10(I.act ^ (<charge>*<charge>) )")

	@Calc:(3,"I.sum = I.sum + ({f<charge>.sum} * (<charge>*<charge>))")
   	@Calc:(3,"chargebalance = chargebalance + ({f<charge>.sum} * (<charge>))")
   
	@Var: f-<charge>.log	0
	@Var: f-<charge>.sum	0
	@Calc:(1,"f-<charge>.sum = 0")
	@Calc:(1,"f-<charge>.log = f<charge>.log")

	@Calc:(3,"I.sum = I.sum + ({f-<charge>.sum} * (<charge>*<charge>))") 
   	@Calc:(3,"chargebalance = chargebalance - ({f-<charge>.sum} * (<charge>))")
}


@Class: davies (i)
This object calculates the activity correction factors for ions charged 1-9 according
to the Davies equation. It also calculates the ionic strength I, and the total sum concentration
of ions with a certain charge.

<br>Hans Meeussen,  June 2002.
{
	@globalvar:  I       <i>     // estimated ionic strength  
	@Var: I.act  1
	@Var: I.sum  0        // sum = sum conc * z^2 for all species
	@Var: I.calc 0        // I as calculated from actual speciation
	@Var: EC     0        // Electrical conductivity millimhos/cm

	@var:  chargebalance 0        // chargebalance of dissolved species
	@Calc:(1, "chargebalance = 0") // initialise chargebalance
	@Calc:(1, "I.act = 10^(0.51* ( (sqrt(I)/(sqrt(I)+1))-(0.23*I)) )"  )


// Modified version according to E. Samson et al Comput. Mat. Science 15 (1999) 285-294
// This modification linearly changes the 0.2 factor from 0.2 at I=0 to 1.5 at I = 1.5.
// which gives a better estimation of ion activity at higher ionic strengths.

//	@Calc:(1, "I.act = 10^(0.51* ( (sqrt(I)/(sqrt(I)+1))-(  (0.2 -  ((I/1.5)*0.05))    * I)) )"  )

      @Calc:(5, "EC      = I.calc/0.013")      // Calculate electrical conductivity from calculated I
	@Calc:(1, "I.sum   = 0") 
	@Calc:(5, "I.sum   = (I.sum -2*I)")
	@Calc:(5, "I.calc  = (I.sum/2)") 

	@Var: f0.log	0
	@Var: f0.sum	0

	@f(1)
	@f(2)
	@f(3) 
	@f(4)
	@f(5)
	@f(6)
	@f(7)
	@f(8)
	@f(9)
} 


@Class: aqphase(){diss}

@Class: species(name, charge){
	@entity(<name>, "@aqphase()", 1)
	@calc_conc_with_davies(<name>, <charge>, "@aqphase()")   
}


@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.
The ".sum" variable is its amount in the min phase.   
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)

November 2004
Hans Meeussen.

{
	@entity(<name>, min) 
	@Var: <name>.un   -1e-3
	@Var: <name>.eq    0
	@Var: <name>.si    0
      @Calc:(2,"<name>.si  = {<name>.logact}/{<name>.si2}") 
	@Calc:(2,"<name>.min = if({<name>.un} >= 0,  ({<name>.un}),   0)  ")
	@Calc:(2,"<name>.eq  = if(0 > {<name>.un},   ({<name>.si} - {<name>.un}), {<name>.si})")
	@Uneq3: unknown:(name:, <name>.un,  delta:, 1e-10, type:, lin, step: , .5) equation:(name:, <name>.eq, tol: , 1e-14, si:, <name>.si )
}



@Class: solid_solution(name, parent_phase)
{
	@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
	// @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
	// 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
//	@Calc:(3, "<name>.eq  = if(0 > <name>.un,   ( log(<name>.<name>) - <name>.un), log(<name>.<name>))")   // log sum fractions
      @Calc:(3, "<name>.eq  = if(0>  <name>.un,   ( log(<name>.<name>) - <name>.un), 1-<name>.<name>)")


	// iterate for a value of the unknown that results in an equation value of zero.

      @Uneq2: unknown:(name:, <name>.un, delta:, 1e-6, type:, lin, step: , .1) equation:(name:, <name>.eq, tol: , 1e-3)
}




@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
   	@Var: 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: cdmusicsurf(name,  parent_phase, surface_area)
{
   @cdmusicsurf(<name>, <parent_phase>, <surface_area>, 1.1, 5)
}


@Class: cdmusicsurf(name,  parent_phase, surface_area, c1, c2)
{
	@globalvar: <name>.c1            <c1>  // default values for capacitance
	@globalvar: <name>.c2            <c2>   // 
	@phase(<name>, <parent_phase>, <surface_area>)
	@phase(<name>_e)
	@link_phase(<name>, <name>_e, 96484.56)

    // The electrostatic layers are child phases of the <name>_e fase, and linked to this phase
    // via the conversion factor  96484.56.
    // This makes the charge balances (<name>_e0.<name>_e.sum) available in coulomb/m2 after phase 4.
	@e_layer (<name>_e0, <name>_e)
	@e_layer (<name>_e1, <name>_e)
	@e_layer (<name>_e2, <name>_e)

    // a diffuse double layer that uses the potential of the e2 layer
	@ddl(<name>, <name>_e2)     

	@Var: <name>.eq1 0.0      // eq1= ddl.sigma + e0.sigma + e1.sigma + e2.sigma
	@Calc:(5,"<name>.eq1 = {<name>.ddl} + {<name>_e0.<name>_e} + {<name>_e1.<name>_e} + {<name>_e2.<name>_e} ")
      @Uneq2: unknown:(name:, <name>_e0.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, <name>.eq1, tol: , 1e-6)


	@Var: <name>.eq2 0.0      // eq2 = (psi0-psi1)* c1 - e0.sigma
	@Calc:(5,"<name>.eq2 = ({<name>.c1}*({<name>_e0.psi}-{<name>_e1.psi}))-{<name>_e0.<name>_e}")
      @Uneq2: unknown:(name:, <name>_e1.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, <name>.eq2, tol: , 1e-6)

	@Var: <name>.eq3 0.0      // eq3 = (psi1-psi2)* c2 - e0.sigma - e1.sigma
	@Calc:(5,"<name>.eq3 = {<name>.c2}*({<name>_e1.psi}-{<name>_e2.psi})-{<name>_e0.<name>_e} - {<name>_e1.<name>_e}")
      @Uneq2: unknown:(name:, <name>_e2.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, <name>.eq3, tol: , 1e-6)
}




@Class: bstern(name, phase, surface_area) 
A Stern model has two electrostatic layers (e0 and e1) 
and a diffuse double layer that uses the potential of the e1 layer.
{
	@Var: <name>.c 3.9 // capacitance in F/m2
	@phase(<name>, <phase>, <surface_area>)

      
	@phase(<name>_e)
	@link_phase(<name>, <name>_e, 96484.56)

	// The electrostatic layers are positioned in the <name>_e fase. This phase has 
	// the sums (<name>_e0.<name>_e.sum)are in the dimension coulomb/m2 after stage 4
	@e_layer (<name>_e0, <name>_e)
	@e_layer (<name>_e1, <name>_e)
	@ddl(<name>, <name>_e1)     

	@Var: <name>.eq1 0.0
	@Calc:(5, "<name>.eq1 = <name>.ddl + <name>_e0.<name>_e.sum + <name>_e1.<name>_e.sum") 
      @Uneq2: unknown:(name:, <name>_e0.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, <name>.eq1, tol: , 1e-6)

	@Var: <name>.eq2 0.0        // eq2 = (psi0-psi1)* c - e0.sigma
	@Calc:(5, "<name>.eq2 = (<name>_e0.psi-<name>_e1.psi) * <name>.c - <name>_e0.<name>_e.sum") 
      @Uneq2: unknown:(name:, <name>_e1.logact, delta:, 1e-6, type:, lin, step: , 1) equation:(name:, <name>.eq2, tol: , 1e-6)

}



//**********************************************************************************
// 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.
This version is adapted to work with the phases objects and now uses name.d.sum rather than 
the old d.name.sum. 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.sum(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 for 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: neutral_diffusion ()
This object class initialises the neutral diffusion calculations, and performs 
part of the calculations that are not component specific.
It is used in combination with the n_diffus() object class in the diffusion.inp file.
<br>Wendy van Beinum, Hans Meeussen
{  @Var: nd.sum1 0     //sum denominator
   @Var: nd.sum2 0     //sum nominator

   //**** Stage 1 ****
   @Calc:(1,"nd.sum1 = 0")    // initialise sums with zero
   @Calc:(1,"nd.sum2 = 0")

   //**** Stage 2 ****
   @Var: nd.factor 0  
   @Calc:(2,"nd.factor = nd.sum1/nd.sum2")
}

@Class: n_diffus (name, charge, D)
This object calculates diffusion for a single component and makes sure that there is no net transport of charge. 
It can be used after a neutral_diffusion() object is defined.
<br>Wendy van Beinum, Hans Meeussen
{  @Var: 1.<name>.diss   0.0
   @Var: 2.<name>.diss   0.0
   @Var: 1.<name>.d  0.0
   @Var: 2.<name>.d  0.0
   @Var:   <name>.z  <charge>
   @Var: 1.<name>.D  <D>

//**** Stage 1 *****
   // dC = C2-C1
   @Var: <name>.dc 0
   @Calc:(1,"<name>.dc = 2.<name>.diss.sum - 1.<name>.diss")

   //Calculate contribution to sum D.z.dC
   @Calc:(1,"nd.sum1 = nd.sum1 + (1.<name>.D * <name>.z * <name>.dc) ")

   //Calculate contribution to sum D.z^2.C
   @Calc:(1,"nd.sum2 = nd.sum2 + (1.<name>.D * <name>.z * <name>.z * 1.<name>.diss)")

//**** Stage 2 ****
   // Calculate J: J = A * D (-dC +z.C. nd.factor)/dx
   @Var: <name>.J  0
   @Var: 1.dx 1
   @Var: 2.A 0
   @Calc:(2,"<name>.J = 2.A * 1.<name>.D*(-<name>.dc + (<name>.z * 1.<name>.diss * nd.factor))/1.dx)")

   //Calculate contribution of J to delta mass of component in both cells

   @Calc:(2, "1.<name>.d = 1.<name>.d - <name>.J")
   @Calc:(2, "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){
     @globalvar: <name>.<phase> <equationvalue>     // define this variable as global (node) variables + default value
     @Uneq2: unknown:(name:, <name>.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, -12) equation:(name:, <name>.<phase>, tol: , 1e-13, default:, <equationvalue>) )
}

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

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

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

// 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: nicamodel(name, parentphase, concentration, donnanVolume){
   // define the nica 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 (act = boltzman factor, sum = 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)    
   @Uneq2: unknown:(name:, <name>_don.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, 0, iia:, true) equation:(name:, <name>_don.<name>, tol: , 1e-4)
   //--------------------------------------------------------------------------------------------
}

@Class: nicamodel(name, parentphase, concentration, donnanVolume, charge){
	@nicamodel(<name>, <parentphase>, <concentration>, <donnanVolume>)
	@Stage:(3,"<name>_don.<name> = {<name>_don.<name>} + <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)

	// The amount of charge that is carried by the empty sites is added to the donnan mass (charge) balance 
	@Stage:(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

	@Stage:(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
	@Stage:(3,"<name>.equation = log10({<name>.<name>} / (10^{<name>.logact})) - log10({<name>.unknown})")

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



@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)
	@calc:(1,"<name>.logact = {f<charge>.log})") 
	@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){
	//@Calc:(1, "<model>_<name>.logact = 0")
	@phase(<model>_<name>, <model>, <density>)
	@entity(<model>_<name>, <model>_<name>, <coef>) 
	@Uneq2: unknown:(name:, <model>_<name>.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, -1, max:, 0) equation:(name:, <model>_<name>.<model>_<name>, tol: , 1e-4, default:, 1)
}

@class: surfsite(model, name, density, coef, charge, plane){
	//@Calc:(1, "<model>_<name>.logact = 0")
	@phase(<model>_<name>, <model>, <density>)
	@entity(<model>_<name>, <model>_<name>, <coef>) 
	@link(<model>_<name>, <model>_<plane>, 0, <charge>)
	@Uneq2: unknown:(name:, <model>_<name>.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, -1, max:, 0) 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: ddl_surface(name, parent_phase, surface_area){
	@phase(<name>, <parent_phase>, <surface_area>)
	@phase(<name>_e)
	@link_phase(<name>,<name>_e, 96484.56)
	@e_layer(<name>_e, <name>_e)
	@ddl(<name>, <name>_e)
	@GlobalVar: <name>.eq 0 // this defines eq as global variables with initial value 0.
	@Stage:(5, "<name>.eq = {<name>.ddl} + {<name>_e.<name>_e}")
      @Uneq2: unknown:(name:, <name>_e.logact,  delta:, 1e-6, type:, lin, step: , 1, default:, 0, iia:, true) equation:(name:, <name>.eq, tol: , 1e-6)
}


@Class: ddl(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_old (name, e_layer){
      // This version of ddl object was used in first Orchestra versions

	@Var: T    298.15          
      @Var: R    8.31441
      @Var: F    96484.56

      // Very old version that uses hard coded expression (Rel[4] type 7)
      //@Var: <name>.ddl 0  
      //@Rel[4]: <name>.ddl <e_layer>.psi 7 1
      //@Calc:(4,"<name>.ddl = I")

      // Modern literal equivalent, with equation in text
	@Var: <name>.ddl 0  
 	@Calc:(4,"<name>.ddl  = if({<e_layer>.psi} > 0, {<name>.ddl}, -{<name>.ddl})")
	@Calc:(4,"<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))   )")    

  
}

@Class: ddl_dzombak (name, e_layer){
      // 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
}



//**
// 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)            
}