#ifndef ISW_HH
#define ISW_HH

// dummy file

#endif // ISW_H

#ifndef SS_HH
#define SS_HH

// dummy file

#endif // SS_H

#ifndef CONSTANTS_HH
#define CONSTANTS_HH

const double	pi	= 3.14159264;
const double	pi2	= pi * pi;
const double	pi4	= pi2 * pi2;

// Reynolds number
const double	ReL	= 2300.0;	// transition from laminar flow
const double	ReT	= 4000.0;	// transition to turbulent flow

#endif // CONSTANTS

#ifndef FADE_HH
#define FADE_HH

#include <math.h>		// tanh

#include "constants.hh"		// pi

inline double
fade(const double x,
     const double x1,
     const double x2,
     const double y1,
     const double y2)
{
  /* fades two functions together smoothly over the range x1 to x2
   * function does not check that x2 > x1
   */
  double theta;
  theta = (x - x1) / (x2 - x1);		// map (linear)     {x1  , x2 } => {0   , +1 }
  theta = (theta - 0.5) * 2.0 * pi;	// map (linear)     {0   , +1 } => {-Pi , +Pi}
  theta = tanh(theta);			// map (non-linear) {-Pi , +Pi} => {-1  , +1 }
  theta = (theta + 1.0) / 2.0;		// map (linear)     {-1  , +1 } => {0   , +1 }

  return (theta * y1 + (1.0 - theta) * y2);
}

inline double
chkfade(const double x,
	const double x1,
	const double x2,
	const double y1,
	const double y2)
{
  double X1 = x1, X2 = x2;
  if (x1 > x2) {
    cerr << "* Warning: chkfade; x2 > x1, swapping" << endl;
    X1 = x2;
    X2 = x1;
  }
  return ((x <= X1) ? y1 : (x > X2) ? y2 : fade(x, X1, X2, y1, y2));
}

#endif // FADE_HH

#ifndef FRICTIONFACTOR_HH
#define FRICTIONFACTOR_HH

#include <iostream>
#include <math.h>
#include <string>

#include "constants.hh"		// ReL, ReT
#include "fade.hh"

inline double
frictionfactor(const double Re, const double r) {
  if (0.0 == Re) {
    return 0.0;
  }
  else if (ReL >= Re) {		// laminar flow
    return 16.0 / Re;		// using k = 4.f.(l/d)
  } else if (ReT <= Re) {	// turbulent flow
    /* S.E.Haaland
     * Simple and explicit formulas for the friction factor in turbulent pipe flow
     * Journal of Fluids Engineering, 105 (1983)
     */
    double A = 6.91 / Re;
    double B = pow((r / 3.71), 1.11);
    double f = pow(-3.6 * log10(A + B), -2);
    return f;
  } else {			// transition region 
    double ffL = frictionfactor(ReL, r);
    double ffT = frictionfactor(ReT, r);
    return fade(Re, ReL, ReT, ffL, ffT);
  }
}

#endif // FRICTIONFACTOR_HH

#ifndef KINEMATICVISCOSITY_HH
#define KINEMATICVISCOSITY_HH

#include <math.h>		// pow
#include <string>

inline double
kerosenekinematicviscosity(const double T)
  /*
   * B.S.Massey
   * Mechanics of fluids
   * ISBN: 0 412 34280 4
   * log-log plot of kinematic viscosity versus temperature is linear for kerosene
   * L(n) = log10(n)
   *
   * T =   0 deg C : nu = 4.0 mm2/s
   * T = 100 deg C : nu = 0.9 mm2/s
   *
   * deg C => K, mm2/s => m2/s
   *
   * T1 = 273.15 : nu1 = 4.0e-6 m2/s
   * T2 = 373.15 : nu2 = 0.9e-6 m2/s
   *
   * L(nu) = m L(T) + c
   *
   *    m = (L(nu2) - L(nu1)) / (L(T2) - L(T1))
   *      = L(nu2/nu1) / L(T2/T1)
   *      = L(0.9/4.0) / L(373.15/273.15)
   *      = -4.781567507
   *
   *    c = L(nu1) - m * L(T1)
   *      = L(4.0e-6) - m * L(273.15)
   *      = 6.251876827 
   *
   * nu {m2/s} = 10^(m * L(T {Kelvin}) + c)
   *
   *    = 10^(m * L(T) + c)
   *    = 10^c * (10^L(T))^m
   *    = 10^c * T^m
   *
   * 10^c = 1.78598097e6
   * 
   * nu = 1.78598097e6 * T^(-4.781567507)
   */
{
  return 1.79e6 * pow(T, -4.78);
}

inline double
kinematicviscosity(const string fluid,
		   const double T)
{
  if ("kerosene" == fluid) {
    return kerosenekinematicviscosity(T);
  } else {
    cerr << __FILE__ << ": fluid \"" << fluid << "\" unknown" << endl;
    exit(-1);
  }
}

#endif // KINEMATICVISCOSITY_HH

#ifndef LIN_HH
#define LIN_HH

#include <iostream>

// translated from lin.cr

// one 2-port, R/C/I; two 2-port, TF/GY
inline double
lin(// parameters
    const causality_t	gain_causality,
    const double	gain,
    // output
    const causality_t	out_causality,
    const int		out_port,
    // input
    const double	input,
    const causality_t	in_causality,
    const int		in_port)
{
  if (out_port == in_port) {			  // R/C/I
    if (gain_causality == in_causality) {
      return input * gain;
    } else {
      return input / gain;
    }
  } else {					  // GY/TF
    if (out_causality == in_causality) {	// gyrator
      if ((out_port == 1 && out_causality != gain_causality)
	  ||(out_port == 2 && out_causality == gain_causality)) {
	return input * gain;
      } else {
	return input / gain;
      }
    } else {				// transformer
      if (out_causality == gain_causality) {
	return input * gain;
      } else {
	return input / gain;
      }
    }
  }
}

// two 2-port, AE/AF
inline double
lin(// parameters
    const double	gain,
    // output
    const causality_t	out_causality,
    const int		out_port,
    // input
    const double	input,
    const causality_t	in_causality,
    const int		in_port)
{
  return
    (out_port == 1) ? input * gain :
    input / gain;
}




// three 2-port, FMR
inline double
lin(// parameters
    const causality_t	gain_causality,
    const double	gain,
    // output
    const causality_t	out_causality,
    const int		out_port,
    // input
    const double	input,
    const causality_t	in_causality,
    const int		in_port,
    const double	modulation,
    const causality_t	mod_causality,
    const int		mod_port)
{
  if (mod_causality == flow) {		// uni-causal
    if (out_port == 2) {
      return 0;
    } else {
      double k = 1.0;
      if (gain_causality == in_causality) {
	k *= gain;
      } else {
	k /= gain;
      }
      if (in_causality == effort) {
	k *= modulation;
      } else {
	k /= modulation;
      }
      return input * k;
    }
  } else {				// bi-causal
    if ((in_causality == effort)
	&&(mod_causality == flow)
	&&(gain_causality == effort)
	&&(in_port == 1)
	&&(mod_port == 1)) {
      return (input / modulation) / gain;
    } else {

	// three 2-port, EMTF
	
	if ((out_causality == gain_causality && out_port == 2)
	    ||(out_causality != gain_causality && out_port == 1)) {
	  return input * gain * modulation;
	} else if((out_causality != gain_causality && out_port == 2)
		  ||(out_causality == gain_causality && out_port == 1)) {
	  return input / (gain * modulation);
	} else {
	  cerr << "* Error: __FILE__ does not cover this case" << endl;
	  exit(-1);
	}
    } // EMTF
  } // bi-causal
}


#endif // LIN_HH

#ifndef PRESSUREDROP_HH
#define PRESSUREDROP_HH

#include <math.h>		// fabs, pow

#include "constants.hh"
#include "frictionfactor.hh"
#include "kinematicviscosity.hh"
#include "sign.hh"

inline double
pressuredrop(const string fluid,
	     const double d,
	     const double l,
	     const double r,
	     const double rho,
	     const double T,
	     const double Q)
{
  double nu = kinematicviscosity(fluid, T);
  double Re = 4.0 * fabs(Q) / (pi * d * nu);
  double f = frictionfactor(Re, r);
  double k = 4.0 * f * l / d;
  double dP = k * 8.0 * rho * pow(Q, 2) / (pi2 * pow(d, 4));
  return (dP * sign(Q));
}

inline double
pressuredrop(const string fluid,
	     const double d,
	     const double l,
	     const double r,
	     const double rho,
	     const double T,
	     const causality_t effort_causality, const int port,
	     const double Q, const causality_t flow_causality, const int port_in)
{
  
  /* assert(effort == causality);
   * assert(flow == causality_in);
   * assert(1 == port_in);
   * assert(1 == port);
   */
  return pressuredrop(fluid, d, l, r, rho, T, Q);
}

#endif // PRESSUREDROP_HH

#ifndef SIGN_HH
#define SIGN_HH

template <class T>
inline int
sign(T x)
{
  return ((x > 0) ? +1 :
	  (x < 0) ? -1 :
	  0);
}

#endif // SIGN_HH

#ifndef SQUARELAW_HH
#define SQUARELAW_HH

#include <math.h>
#include "sign.hh"

inline double squarelaw(const double gain,
			const causality_t causality, const int port,
			const double input, const causality_t in_causality, const int in_port)
  /*
   * implements P = R Q^2
   * direction of flow is retained
   */
{
  if (causality == effort) {
    return pow(input, 2) * gain * sign(input * gain);
  } else {
    return sqrt(fabs(input / gain)) * sign(input / gain);
  }
}

#endif // SQUARELAE_HH

#ifndef STATICPRESSURE_HH
#define STATICPRESSURE_HH

#include <math.h>		// log, pow

#include "constants.hh"

inline double
staticpressure(const double beta,
	       const double C_d,
	       const double d,
	       const double P_ref,
	       const double rho,
	       const causality_t causality, const int port,
	       const double Q1, const causality_t causality1, const int port1,
	       const double Q2, const causality_t causality2, const int port2)
{
  static double P;
  if (0.0 != Q1 && 0.0 != Q2) {
    double num = pi2 * pow(d, 4) * log(Q1 / Q2);
    double den = 8.0 * beta * rho * Q1 * (Q2 - Q1 + C_d * (Q1 + Q2) / 2.0);
    P = P_ref + log(num / den)/beta;
  }
    return P;
}

#endif // STATICPRESSUE_HH

%SUMMARY CT2    Constitutive Relationship for a two port thermo C
%DESCRIPTION Parameter 1: c_v (specific heat at constant volume)
%DESCRIPTION Parameter 2: gamma = c_p/c_v
%DESCRIPTION Parameter 3: mass of (ideal) gas within component.
%DESCRIPTION Parameter 4: t_0 -- the temperature at which internal
%DESCRIPTION energy is zero.


%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  1997/12/07 20:45:21  peterg
% %% Initial revision
% %%
% %% Revision 1.1  1996/11/02  10:21:19  peterg
% %% Initial revision
% %%
% %% Revision 1.1  1996/09/12 11:18:26  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%SUMMARY StefanBoltzmann: Stefan-Boltzmann radiation law.	
%DESCRIPTION Parameter 1: Stefan-Boltzmann constant


%DESCRIPTION Parameter 2: Area of radiating surface



















OPERATOR StefanBoltzmann;










FOR ALL sigma,Area,input


LET StefanBoltzmann(sigma,Area,flow, 1, 

	input, effort, 1)









	 = sigma*area*input^4;

%SUMMARY cm: relation for 2-port CM component
%DESCRIPTION Parameter 1 capacitance at separation x_0
%DESCRIPTION Parameter 2 x_0
%DESCRIPTION parameter 3 moving-plate mass

%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  1996/11/02  10:21:19  peterg
% %% Initial revision
% %%
% %% Revision 1.1  1996/09/12 11:18:26  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


OPERATOR cm;
%Linear electrical bit
FOR ALL c_0,x_0,elec_state,mech_state LET
cm(c_0,x_0,effort,1,
	elec_state,state,1,
	mech_state,state,2
	)
	= elec_state/(c_0*x_0/mech_state);

%Nonlinear mechanical bit
FOR ALL c_0,x_0,elec_state,mech_state LET
cm(c_0,x_0,effort,2,
	elec_state,state,1,
	mech_state,state,2
	)
	= -(c_0*x_0)*((elec_state/mech_state)^2)/2; 

END;;

%SUMMARY cr generic CR
%DESCRIPTION Argument is an algebraic expression with no embeddedwhite space
%DESCRIPTION Only available for one ports just now
%DESCRIPTION effort (or integrated effort) variable must be called mtt_e
%DESCRIPTION flow (or integrated flow) variable must be called mtt_f
%DESCRIPTION For example:
%DESCRIPTION             mtt_e=k*mtt_f
%DESCRIPTION             mtt_f=mtt_e/r

% $Log$
% Revision 1.3  2000/10/05 10:13:00  peterg
% New eqn2ass function.
% Started extension to multiports
%
% Revision 1.2  2000/10/03 18:35:04  peterg
% Removed comment bug
%
% Revision 1.1  2000/10/03 18:34:00  peterg
% Initial revision
%

% Flow output
FOR ALL mtt_cr_e,mtt_cr_f, input, in_cause
LET cr(mtt_cr_e,mtt_cr_f,flow, 1, input, in_cause, 1) 
    = sub(mtt_e=input,mtt_cr_e);


%%% Q&D FMR 2 port.
FOR ALL mtt_cr_e,mtt_cr_f,input_1,input_2
LET cr(mtt_cr_e,mtt_cr_f,flow,1,
	input_1,effort,1,
	input_2,flow,2
	)  = sub(mtt_mod=input_2,sub(mtt_e=input_1,mtt_cr_e));

END;

%SUMMARY delta_h	CR for gas turbine compressor


OPERATOR delta_h;

% Port 1 - generates delta h
FOR ALL c_p,Temperature,Massflow,DeltaT
LET delta_h(c_p, flow, 1,
		Temperature,effort,1,
		Massflow,flow,2,
		DeltaT,effort,3)
	 = Massflow*c_p*DeltaT;

% Port 2 - generates zero effort
FOR ALL c_p,Temperature,Massflow,DeltaT
LET delta_h(c_p, effort, 2,
		Temperature,effort,1,
		Massflow,flow,2,
		DeltaT,effort,3)
	 = 0;

% Port 3 - generates zero effort
FOR ALL c_p,Temperature,Massflow,DeltaT
LET delta_h(c_p, flow,3,
		Temperature,effort,1,
		Massflow,flow,2,
		DeltaT,effort,3)
	 = 0;

%SUMMARY lin    linear constitutive relationship
%DESCRIPTION Parameter 1 defines input causality relating to parameter 2
%DESCRIPTION value is effort, flow or state
%DESCRIPTION Parameter 2 is the gain corresponding to the causality of
%DESCRIPTION parameter 1.
%DESCRIPTION Supported components:

%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linear constitutive relationship.


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % $Id$
% % $Log$
% % Revision 1.3  1998/07/04 10:47:04  peterg
% % back under RCS
% %
% % Revision 1.2  1998/03/04 15:38:54  peterg
% % Added END statement
% %
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%DESCRIPTION    single port components: R,C,I
%Linear Constitutive Relationship for single port components: R,C,I.
% e = Gain*f (if gain_causality = flow) 
%           f = Gain*e (if gain_causality = effort)
OPERATOR lin;
FOR ALL gain_causality, gain, causality, input, other_causality
SUCH THAT causality = gain_causality
LET lin(gain_causality, gain, other_causality, 1, input, causality, 1)
         = gain*input;

%Linear CR: e = (1/Gain)*f (if gain_causality = flow) 
%           f = (1/Gain)*e (if gain_causality = effort)
FOR ALL gain_causality, gain, causality, input, other_causality
SUCH THAT causality NEQ gain_causality
LET lin(gain_causality, gain, other_causality, 1, input, causality, 1)
         = input/gain;

%DESCRIPTION    two port components: AE, AF
% Linear Constitutive Relationship for AE and AF
% Output = gain * input

% Unicausal form
FOR ALL gain, input, causality
LET lin(gain, causality, 2, input, causality, 1) = gain*input;

%Bicausal form
FOR ALL gain, output, causality
LET lin(gain, causality, 1, output, causality, 2) = output/gain;

%DESCRIPTION   two port component: TF
% Linear Constitutive Relationship for TF
FOR ALL gain_causality, gain, causality, outport, input, same_causality, inport

SUCH THAT 
       ( causality = same_causality ) 
       AND 
       ( inport NEQ outport )
       AND
       (
       ( (causality = gain_causality) AND (outport = 2) )
       OR
       ( (causality NEQ gain_causality) AND (outport = 1) )
       )
LET lin(gain_causality, gain, causality, outport, 
          input, same_causality, inport)
        = gain*input;

FOR ALL gain_causality, gain, causality, outport, 
        input, same_causality, inport
SUCH THAT 
       ( causality = same_causality ) 
       AND       
       ( inport NEQ outport )
       AND
       (
       ( (causality NEQ gain_causality) AND (outport = 2) )
       OR
       ( (causality = gain_causality) AND (outport = 1) )
       )
LET lin(gain_causality, gain, causality, outport, 
         input, same_causality, inport)
        = input/gain;

%% This version in not reliable. I rellly need to pass component names
%% as cr arguments.

%DESCRIPTION    two port component: GY
% Linear Constitutive Relationship for GY

FOR ALL gain, input, causality, gain_causality, other_causality, 
        outport, inport
SUCH THAT 
        (causality NEQ other_causality) 
        AND
        ( inport NEQ outport )
        AND
        (
        ( (causality NEQ gain_causality) AND (outport = 2) )
        OR
        ( (causality NEQ gain_causality) AND (outport = 1) )
        )
LET lin(gain_causality, gain, other_causality, outport, 
        input, causality, inport)
         = input/gain;

FOR ALL gain, input, causality, gain_causality, other_causality, 
        outport, inport
SUCH THAT 
        (causality NEQ other_causality) 
        AND
        ( inport NEQ outport )
        AND
        (
        ( (causality = gain_causality) AND (outport = 2) )
        OR
        ( (causality = gain_causality) AND (outport = 1) )
        )
LET lin(gain_causality, gain, other_causality, outport, 
        input, causality, inport)
         = gain*input;

%DESCRIPTION    three port component: FMR

% Linear Constitutive Relationship for FMR - unicausal case
% Flow modulation multiplies effort on port 1 (or divides flow)

% The 4 possibilities follow...
FOR ALL gain_causality, gain, out_causality, input, in_causality,
        mod_input
SUCH THAT (gain_causality=in_causality) AND (out_causality=flow)
LET lin(gain_causality, gain, out_causality, 1, 
                input, in_causality, 1,
                mod_input, flow, 2)
         = input*gain*mod_input;

FOR ALL gain_causality, gain, out_causality, input, in_causality,
        mod_input
SUCH THAT (gain_causality=in_causality) AND (out_causality=effort)
LET lin(gain_causality, gain, out_causality, 1, 
                input, in_causality, 1,
                mod_input, flow, 2)
         = input*gain/mod_input;

FOR ALL gain_causality, gain, out_causality, input, in_causality,
        mod_input
SUCH THAT (gain_causality NEQ in_causality) AND (out_causality=flow)
LET lin(gain_causality, gain, out_causality, 1, 
                input, in_causality, 1,
                mod_input, flow, 2)
         = input*mod_input/gain;

FOR ALL gain_causality, gain, out_causality, input, in_causality,
        mod_input
SUCH THAT (gain_causality NEQ in_causality) AND (out_causality=effort)
LET lin(gain_causality, gain, out_causality, 1, 
                input, in_causality, 1,
                mod_input, flow, 2)
         = input/(gain*mod_input);

% Linear Constitutive Relationship for FMR - bicausal case
% Deduces the flow on port 2.

% The 2 possibilities follow...
FOR ALL gain,  e_input, f_input
LET lin(effort, gain, flow, 2, 
                e_input, effort, 1,
                f_input, flow, 1)
         = (f_input/e_input)/gain;

%EMTF component - modulation only
% Linear Constitutive Relationship for EMTF
FOR ALL gain_causality, gain, causality, outport, input, same_causality, inport
SUCH THAT 
       ( (causality = gain_causality) AND (outport = 2) )
       OR
       ( (causality NEQ gain_causality) AND (outport = 1) )
LET lin(gain_causality, causality, outport, 
          input, same_causality, inport,
	  gain, effort, 3)
        = gain*input;

FOR ALL gain_causality, gain, causality, outport, 
        input, same_causality, inport
SUCH THAT 
       ( (causality NEQ gain_causality) AND (outport = 2) )
       OR
       ( (causality = gain_causality) AND (outport = 1) )
LET lin(gain_causality, causality, outport, 
         input, same_causality, inport,
	  gain, effort, 3)
        = input/gain;

%EMTF component - modulation and gain
% Linear Constitutive Relationship for EMTF
FOR ALL gain_causality, gain, causality, outport, input,
same_causality, inport, modulation
SUCH THAT 
       ( (causality = gain_causality) AND (outport = 2) )
       OR
       ( (causality NEQ gain_causality) AND (outport = 1) )
LET lin(gain_causality, gain, causality, outport, 
          input, same_causality, inport,
	  modulation, effort, 3)
        = gain*modulation*input;

FOR ALL gain_causality, gain, causality, outport, 
        input, same_causality, inport, modulation
SUCH THAT 
       ( (causality NEQ gain_causality) AND (outport = 2) )
       OR
       ( (causality = gain_causality) AND (outport = 1) )
LET lin(gain_causality, gain, causality, outport, 
         input, same_causality, inport,
	  modulation, effort, 3)
        = input/(gain*modulation);

END;;

%%% linx - cr for single port I and C with an initial state x_0



%DESCRIPTION    linear cr for single port I and C with an initial state x0
%DESCRIPTION    only adds x0 if in integral causality


OPERATOR linx;

%% Input causality as specified
%Linear Constitutive Relationship for single port components: C,I.
% e = Gain*f (if gain_causality = flow) 
%           f = Gain*e (if gain_causality = effort)

FOR ALL gain_causality, gain, causality, input, other_causality
SUCH THAT (causality = gain_causality) AND (causality = state)
LET linx(gain_causality, gain, x0, other_causality, 1, input, causality, 1)
         = gain*(input + x0);


FOR ALL gain_causality, gain, causality, input, other_causality
SUCH THAT (causality = gain_causality) AND (causality NEQ state)
LET linx(gain_causality, gain, x0, other_causality, 1, input, causality, 1)
         = gain*(input);

%% Input causality not as specified
%Linear CR: e = (1/Gain)*f (if gain_causality = flow) 
%           f = (1/Gain)*e (if gain_causality = effort)

FOR ALL gain_causality, gain, x0, causality, input, other_causality
SUCH THAT (causality NEQ gain_causality) AND  (causality = state)
LET linx(gain_causality, gain, x0, other_causality, 1, input, causality, 1)
         = (input+x0)/gain;

FOR ALL gain_causality, gain, x0, causality, input, other_causality
SUCH THAT (causality NEQ gain_causality) AND  (causality NEQ state)
LET linx(gain_causality, gain, x0, other_causality, 1, input, causality, 1)
         = (input)/gain;




END;;

%SUMMARY oneway	One way constitutive relationship eg Diode
%DESCRIPTION Parameter 1 is a large number being the forward gain
%DESCRIPTION	-- the reciprocal is the backward gain
%DESCRIPTION The input must be an effort
%DESCRIPTION Typical use is an R component with effort input

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


OPERATOR oneway; 

%Input has flow causality
FOR ALL r,  input
LET oneway(r, effort, 1, input, flow, 1)
	 = ((1 - sign(input))/2)*r*input;

%Input has effort causality
FOR ALL r,  input
LET oneway(r, flow, 1, input, effort, 1)
	 = ((1 - sign(input))/2)*(1/r)*input;

%SUMMARY polytrop	CR for gas turbine compressor


OPERATOR polytrop;

% Port 1 generates zero flow
FOR ALL Ipressure,temperature,Fpressure,gamma,enthflow
LET polytrop(gamma, flow, 1,
		Fpressure,effort,1,
		enthflow,flow,2,
		temperature,effort,3,
		Ipressure,effort,4)
	 = 0;

% Port 2 generates deltaT
FOR ALL Ipressure,temperature,Fpressure,gamma,enthflow
LET polytrop(gamma, effort, 2,
		Fpressure,effort,1,
		enthflow,flow,2,
		temperature,effort,3,
		Ipressure,effort,4)
	 = temperature*((Ipressure/Fpressure)^(gamma)-1);

% Port 3 generates zero flow
FOR ALL Ipressure,temperature,Fpressure,gamma,enthflow
LET polytrop(gamma, flow, 3,
		Fpressure,effort,1,
		enthflow,flow,2,
		temperature,effort,3,
		Ipressure,effort,4)
	 = 0;

% Port 4 generates zero flow
FOR ALL Ipressure,temperature,Fpressure,gamma,enthflow
LET polytrop(gamma, flow, 4,
		Fpressure,effort,1,
		enthflow,flow,2,
		temperature,effort,3,
		Ipressure,effort,4)
	 = 0;

%SUMMARY powerlaw	powerlaw constitutive relationship
%DESCRIPTION Parameter 1 defines input causality relating to parameter 2
%DESCRIPTION value is effort or flow
%DESCRIPTION Parameter 2 is the gain r corresponding to the causality of
%DESCRIPTION parameter 1.
%DESCRIPTION Supported components:


%DESCRIPTION	single port components: R

%Powerlaw Constitutive Relationship for single port components: R


OPERATOR powerlaw;
FOR ALL gain_causality, gain, power, causality, input, other_causality
SUCH THAT causality = gain_causality
LET powerlaw(gain_causality, gain, power, other_causality, 1, input, causality, 1)
	 = gain*(abs(input)^power)*sign(input);


FOR ALL gain_causality, gain, power, causality, input, other_causality
SUCH THAT causality NEQ gain_causality
LET powerlaw(gain_causality, gain, power, other_causality, 1, input, causality, 1)
	 = ( (abs(input)/gain)^(1/power) )*sign(input);

END;

%SUMMARY reed Nonlinear 2-port R for musical reed component


%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%     %%%%% Model Transformation Tools %%%%%
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linear constitutive relationship.

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %% Revision 1.1  1996/11/02 10:21:19  peterg
% %% Initial revision
% %%
% %% Revision 1.1  1996/09/12 11:18:26  peter
% %% Initial revision
% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Linear Constitutive Relationship for reed - unicausal case with 
% pressure output.
OPERATOR reed,abs,sign;

% Port 1 - the modulated R
FOR ALL D,q,H,airflow,displacement
LET reed(D,q,H, effort, 1, 
		airflow,flow,1,
		displacement, effort,2)
	 = (D*sign(airflow)*(airflow)^q)/((H-displacement)^2);

% Port 2 - zero flow
FOR ALL D,q,H,airflow,displacement
LET reed(D,q, flow, 2, 
		airflow,flow,1,
		displacement, effort,2)
	 =0;

%SUMMARY sat_tank Saturation nonlinearity with variable slopes for tank
%DESCRIPTION Parameter 1 is the slope of the "normal" linear part of the CR
%DESCRIPTION Parameter 2 is the "large" slope
%DESCRIPTION Parameter 3 is the lower bound of the state
%DESCRIPTION Parameter 4 is the upper bound of the state


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Version control history
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% $Id$
% %% $Log$
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

OPERATOR sign0;
FOR ALL x
	LET sign0(x) = (sign(x)+1)/2;

OPERATOR sat_tank; 
%Output has effort causality, input is state
FOR ALL k_0, k_1, x_0, x_1, x
LET sat_tank(k_0, k_1, x_0, x_1, effort, 1,
         x, state, 1)
	 = x*k_0 + (x-x_1)*(k_1-k_0)*sign0(x-x_1) + (x-x_0)*(k_1-k_0)*sign0(x_0-x);

%DESCRIPTION Sensitivity version of lin

OPERATOR slin;
FOR ALL gain_causality, gain, causality, input, other_causality
LET slin(gain_causality, gain, other_causality, 1, input, causality, 1)
         = lin(gain_causality, gain, other_causality, 1, input,
	 causality, 1);
END;;

%SUMMARY square	square-law constitutive relationship
%DESCRIPTION Parameter 1 defines input causality relating to parameter 2
%DESCRIPTION value is effort or flow
%DESCRIPTION Parameter 2 is the gain r corresponding to the causality of
%DESCRIPTION parameter 1.
%DESCRIPTION Supported components:


%DESCRIPTION	single port components: R

%Square-Law Constitutive Relationship for single port components: R
% output = Gain*input^2*sign(input) {if gain_causality = causality} 
% output = (1/Gain^(1/2))*input^(1/2)*sign(input) 
%		{if gain_causality not= causality} 



OPERATOR square;
FOR ALL gain_causality, gain, causality, input, other_causality
SUCH THAT causality = gain_causality
LET square(gain_causality, gain, other_causality, 1, input, causality, 1)
	 = gain*input^2*sign(input);


FOR ALL gain_causality, gain, causality, input, other_causality
SUCH THAT causality NEQ gain_causality
LET square(gain_causality, gain, other_causality, 1, input, causality, 1)
	 = input^(1/2)*sign(input)/(gain^(1/2));


END;