TITLE Example of model with definite integral
: p(q) contains a definite integral plus
:   an explicit algebraic term.
: The model also contains two examples of ratios
:   that can go to 0/0 .
: Default output - the variable q decreases monotonically
:   from 1 in three different curves, for L = .2, 1, and 5
:
: Simulation Resources, Inc.
: M. Kootsey
: 3/7/2001

STEPPED
{
  L = 0.2,1.0,5.0
}

INDEPENDENT
{
  q FROM 0 TO 10 WITH 100
}

ASSIGNED
{
  p FROM 0 TO 1.1
  z
  delta_q
}

STATE
{
  y
}

FUNCTION sinratio(r)
{                : We have to provide for the limit when r is near zero
  IF(r < 1.e-4)
  {
    sinratio = 1.0
  } ELSE
  {
    sinratio = sin(r) / r
  }
}

DERIVATIVE term
{
  y' = sinratio(z)  : This is the integrand in the definite integral
}

FUNCTION integral(qq)  : This function does the definite integral 
{                      : using the differential equation solver
LOCAL saveq

  IF(qq > 1.e-4)
  {
    saveq = q  : We have to save the main independent variable here
    q = 0.0    : and restore it below. q becomes a temporary
               : independent variable for the definite integral.
    delta_q = qq * L / 100. : qq sets the limit for the integral
    y = 0.0    : Initialize the variable to be integrated
    FROM j = 0 TO 100 BY 1    : This loop does the definite integral
    {
      z = j * qq * L / 100.
      SOLVE term
    }
    integral = 2. * y /(qq * L)
    q = saveq  : Restore the main independent variable
  } ELSE
  {
    integral = 2. / L : This is the limit in the 0/0 case
  }
}

PLOT p VS q

BREAKPOINT
{
  p = integral(q) - pow(sinratio(q * L/2.), 2.)
}