CONTROL help

• Introduction
• DMI
  • DMI SET
  • DMI PRESET, DMI TIME
  • DMI GO, DMI ON, DMI HLT
  • DMI : test and service tools
  • DMI : CIO variables created and used
• MTR
  • MTR INQuire
  • MTR INIT
  • MTR POS, MTR RATE
  • MTR Goto
  • MTR Check
  • MTR BRake
  • MTR DISable
  • MTR TEST
  • MTR : errors and warnings
• XYPOS
• NMR
  • MTR CAL
  • MTR ION
  • MTR E, NMR F
• EXPORT
• COCON
• Goniometry
  • GON INIT
  • GON MS
  • GON SM
  • SECT
  • DIPS
  • CHAN2
  • CHAN3
• Reaction kinematics
  • KIN
  • DEPth
  • ISOtopes
  • MASS
• Spectrum simulation
  • RBSIM
• Fit, smooth, interpolate
  • LLSQ
  • SPLINE
  • SPLINT

Introduction

    Command line parameters:
      -h                display list of commands
      -s spectrumfile   (LU 1, read only)
      -a spectrumfile   (LU 2, read/write)
      -d asciifile      (LU 7)
      -l logfile        (LU 14)
      -p cio_program    (LU 15, for command RUN)

    Each of this options should be given separately, each starting
    with " -".
    Further command line arguments are treated as normal commands.

DMI

    dmi set   : set all sorts of parameters
    dmi preset: preset time or number of clock pulses
    dmi time  : print out time and clock pulses of current measurement
    dmi go    : start a measurement, go in wait loop while measuring.
    dmi on    : start a measurement, no builtin wait loop
    dmi hlt   : stop a measuremnt
    dmi load/debug/flush/regs : debugging tools

DMI SET

    dmi set <number of ADC's> <ADC numbers> <OPTION> <additional input>

    OPTION :
      DMI           collect singles spectra
      DMICL         collect singles spectra, start on first real time
                    clock pulse.
      COIN,COINCL   coincident mode, for XY position sensitive detector
      DMA,DMACL     list mode.
      MULTIscale    multiscaling mode. ADC's work as for DMI, but a
                    spectrum of yield versus time is also generated.
    additional input varies for different options

DMI PRESET, DMI TIME

    dmi preset : preset time or number of clock pulses
    dmi time   : print out time and clock pulses of current measurement

    clock time (specified in seconds)
    real 'time' = number of pulses fed to multiplexer interface input
    live 'time' (one counter for each selected ADC)
                = number of 'real' time pulses, but filtered through
                  a gate that closes when the ADC is busy

DMI GO, DMI ON, DMI HLT

    dmi go  : start a measurement, go in wait loop while measuring.
    dmi on  : start a measurement, no builtin wait loop
    dmi hlt : stop a measuremnt

DMI : test and service tools

    dmi load  : down-load adc multiplexer
    dmi flush : flush output from adcmultiplexer to OS-90
    dmi regs  : print adc multiplexer internal registers
    dmi debug : control output to logfile /dd/tmp/dmidma.log
                asks for an integer number DEBUG
                if (DEBUG == 0) no logfile output
                if (DEBUG & 1) very limited output
                if (DEBUG & 2) report every buffer switch
                if (DEBUG & 4) report every coincident event

DMI : CIO variables created and used

    int ADCBUSY         0 when stopped, 1 when measuring.
    double dmierror     counts number of errors detected by the program
    int prekloks[9]     preset number of clock pulses
    int adckloks[9]     elapsed - preset number of clock pulses
    double adctime[8]   elapsed times, set by DMI TIME or end of
                        measurement
    double dmipid       set by command DMI CHECK
                        <= 0 if control is offline (DMI not permitted)
    int adcabc[]        data used in COIN or MULTIscale mode
                        in multiscaling mode:
                        adcabc[0] : current channel number in time
                                    spectrum
                        adcabc[1] : time interval in clock ticks
                                    (units of 1/256 s)

    0: number of buffer switches
    1: number of events
    2: number of incomplete events
    3: number of other syntax errors in events
    4: number of events with x or y out of range
    5: number of events with x y and z equal to zero

MTR

MTR INQuire

    MTR INQuire <devicename>
    
      This should always be the first MTR command of a session.
      When the device name is known to the program already,
      this command may be abbreviated to

        MTR INQ;

      Device name is the OS9-system RS232 port (/t1). The current
      status and position of the ISMC are reported. Therefore this
      is also a useful command after any errors (in addition to
      MTR Check, see below). Note that MTR INQ also works when the
      motors are not initialized, while MTR Check assumes that
      they have been initialized.

MTR INIT

    MTR INIT <parameters>

      This command asks for the following parameters:
        Number of motors (2 or 3}
        Number of encoder bits (13 or 14) for each motor
        Number of motor steps / 360 degrees
        Number of encoder increments / 360 degrees
        Encoder value where the angle is 0
        Minimum and maximum angle (allowed range)
        Undershoot and required precision (both in degrees)
        Current angles (degrees, may also be set by MTR POS)
        Parameters R,S,F,Z, determining the rate and the acceleration
          algorithm (may also be set by MTR Rate)
          R is top speed (for R < 41 the rate is 100*R steps/second)
          F is first rate, i.e. the minimum speed.
          S and Z determine the acceleration from "F" to "R" or vice versa
          the acceleration is proportional to 1/{Z*(256-S)}.

MTR POS, MTR RATE

MTR Goto

    MTR Goto <angles>

      Arguments of this command are Nmotor angles (in degrees).

MTR Check

    MTR Check

      Besides the angle values based on the absolute encoders, the
      OS90-system and the ISMC each keep track of what the angles
      should be based on the the number of motor steps. This
      'ISMC angle' is only initialized by MTR INIT and MTR POS.
      The OS90-system angle is adjusted to the encoder value
      before each step request (ANG GO or ANG SCAN). The 
      (hopefully not too) different values are reported by
      MTR CHECK.

      MTR CHECK will also report on current status errors.

MTR BRake

    MTR BRake

      May be useful to abort an erroneous ANG GO or ANG SCAN command.
      Stops all motors. A new MTR POS or MTR INIT is then required.

MTR DISable

    MTR DISable

      This command should be given at the end of a session, and
      makes the motors inaccessible until a new command MTR INIT
      is given.

MTR TEST

    MTR TEST <keyword> [<text>]

      Possible keywords are:
        TRON : trace all input/output to the ISMC
        TROF : turn tracing off
        WR text: write text to the ISMC, and read response

MTR : errors and warnings

    MTR POS;

XYPOS

    XYPOS x, y

    Move the XY-table to position (x,y). Default is current position.

NMR

    A1=0.1747515
    A2=0.3195214D-9

NMR CAL

    NMR CAL <isotope1> <energy1> <frequency1> <isotope2> <energy2> <frequency2>

      Calculates the coefficients a1 and a2 from two calibration points.
      Isotopes specified as mnemonics, e.g. P (=proton), ALPHA, 4HE.
      Energies in keV, frequencies in kHz. To obtain the current calibration
      type

        NMR CAL P 0 0

      To obtain a linear calibration, set fequency2 = 0.

NMR ION

    NMR ION <charge state> <Nisotopes> <Nisotopes*(Isotope, Ni)>


      charge state : the charge state of the ion
      Nisotopes :    the number of different isotopes in the ion
      Isotope :      mnemonic, e.g. P (=proton), 14N
      Ni :           number of atoms of this isotope per ion

      This command specifies the ion beam.
      Example: to specify a beam of HHD+ ions give
      
        NMR ION 1,2, P,2, D,1

      The command NMR ION has to be given only once after loading CONTROL, or
      whenever the beam ion is changed. There is no default ion.

NMR E, NMR F

    NMR E <energy>
    NMR F <frequency>

      energy :    in keV
      frequency : in kHz
      These commands relate the energy of the (first) nucleus of the ion to
      the NMR frequency. Results are listed on the console and also assigned
      to variables R0=frequency and R1=energy.
      They do not affect the NMR control of the Van de Graaff.

EXPORT

    The command EXPORT should only be given after a session has been initialized
    and all relevant arrays have been defined. The command EXPORT enables the
    current set of CIO variables to be shared with other copies of program
    control running in the same computer. As a side effect, all current CIO
    variables effectively become permanent. They can not be deleted, although
    their value may still be modified.

    This command is only available in the OS-9 version of control.

    Note:
    To enable deletion of variables after giving export, proceed as follows.
      a) stop 'slave' programs, started by COCON.
      b) give the command 'keep -z'
      c) delete the relevant variables
         ...
      d) give the command export
      e) Now you may restart cocon.

COCON

    The OS-9 command COCON starts up a version of program CONTROL that 'grabs'
    the CIO variables exported (see command EXPORT) by another copy of control.
    In one OS-9 computer several copies of COCON may be started simultaneously.
    Such a program may be used to display spectra, inspect data, and even
    modify some global parameters, such as those set by DMI PRESET.

    This command is only available in the OS-9 version of control.

Goniometry

    0 <= theta <= 180, -180 <= phi <= 180.

GON INIT

    GON Init <theta0> <phi0> <theta1> <phi1>

      This call specifies frame S with respect to frame M.
        theta0, phi0 : coordinates of the z' axis (theta' = 0)
        theta1, phi1 : coordinates of some point of the (z',+x') plane (phi'=0)

    The parameters given via GON Init are stored in a common block and may be
    used for many further calls GON MS and GON SM. GON MS and GON SM store
    their results in CIO variables
      Gonphi   : theta' or theta
      Gontheta : phi' or phi

GON MS

    GON MS <theta> <phi>

      Calculates theta',phi'

GON SM

    GON SM <theta'> <phi'>

      Calculates theta,phi

SECT

    SECT <theta1> <phi1> <theta2> <phi2> <theta3> <phi3> <theta4> <phi4>

      The input angles specify four points on a unit sphere.
      Let points 1 and 3 define one great circle, and points 2 and 4 another one.

      SECT calculates the point of intersection of the two circles and also the
      angle between the two planes. Of the 2 possible solutions, the one closest
      to point 1 is returned. All input and output is in degrees.

DIPS

    Uses the output of CIO program CIRCLE, TRIANGLE or SQUARE.
    The program asks a string (Circle, Square or Triangle), the begin point of
    the circle, square or triangle. and a dip recognition criterion.
    More dips are accepted when this number is higher.(normal values: 0.5-2).

    The circle, square or triangle is drawn on the screen with the dips.
    The dips are numbered in order of magnitude. Then the numbers are asked
    of dips one expects to belong to the same plane. The plane is drawn and the
    intersection with the other planes is calculated

CHAN2

    This is a program to do goniometry for the 2-axis RBS goniometer.
    The two motor angles are called (theta,phi) in the program.

    CHAN calculates the directions (theta,phi) of all strings and planes with
    specified indices. For calibration two strings or one string and a plane
    have to be given first. These calibration strings or string + plane are to
    be input as H,K,L, theta,phi (in degrees).

    If two strings are given the order and sign of the indices should be such
    that the angle between the strings is consistent with the given directions.
    If the difference is more then 4 degrees an error message is given and new
    input is asked. If the error is less , the program can perform a correction
    by assuming an offset in theta, i.e. theta is replaced by (theta - th0) in
    the calculation. This only makes sense if this offset is the cause of the
    deviation

    The program checks if there is a unique solution. If there are two solutions
    both are calculated.

CHAN3

    This is a program to do goniometry for the 3-axis LEIS goniometer.
    For every set of motor angles there is a certain orientation of the crystal.
    To calibrate the relation between these, the program asks for two axial
    directions, (or for one axial direction and a point of a plane,) and
    corresponding motor angles. These are sets of angles where the specified
    axial direction (or plane) coincides with the beam direction. The axial
    directions and planes are specified by their Miller indices h,k,l. Once
    this relation is established, there are several options:

    Order :
      Input a positive integer, and a value for motor angle 3.
      The program will calculate and print all accessible strings with one index
      equal to 'order', and for each of these the motor angles, and the angle
      of the normal of the crystal with up to 3 detectors, and with the beam.

    Hkl :
      Input the required orientation of the crystal, by specifying two directions,
      and corresponding angles. Calculate quantities of interest.

    Motpos :
      Input motor angles. Calculate quantities of interest.

    Rotate :
      Input initial motor angles, an axis (h,k,l) of the crystal, and the angle
      of rotation. Calculate the new motor positions to effectuate the requested
      rotation.

    The lab frame used in the program has its X axis vertical, Y axis
    perpendicular to the beam, and Z axis along the beam. Polar coordinates
    theta (=angle with Z axis) and phi (=angle with ZX plane) are also used.
    Thus the horizontal plane has phi= +90 or -90 degrees.

    In subsequent calls of CHAN3 one can skip the calibration section and go
    directly to the options section by typing:

      chan3;

Reaction kinematics

KIN

    Relativistic reaction kinematics of 2-particle reactions.

                m3
               /
              / psi
     m1----->m2----
              \ zeta
               \
                m4

    The reaction is m1 + m2 --> m3 + m4, where
      m1 : is beam particle (may also be GAMMA)
      m2 : is target particle, originally at rest
      m3 and m4 : are reaction products

    OPTION INPUT             OUTPUT
      1    1 value of psi    T3,T4,V3,V4,psi,theta,zeta,cmlab
      2    2 values of psi   same quantities averaged
      3    N values of psi   specified quantities for N angles
    4,5,6  same as 1,2,3 but theta values are input
      0    STOP

    One may either specify the lab. angle psi (options 1 and 2), or the
    center-of-mass angle theta (options 3 and 4), of particle m3.
    Several output options are available. Output quantities may be
      kinetic energies : T1,T2,T3,T4
      momenta (MeV/c)  : P1,P2,P3,P4
      velocities (v/c) : V1,V2,V3,V4
      lab angles PSI (for m3) and ZETA (for m4), c.m. angle THETA,
      ratio of cross sections cm/lab for particles m3
    Particles m1,m2,m3,m4 are specified by their mnemonics, such as
    27AL, 4HE (or ALPHA, or A), GAMMA (or G).
    The requested input is explained when the program is run.

DEPth

    Calculates energy-depth (or depth-channel) relations, for particles
    leaving a material after elastic scattering or after a nuclear reaction,
    or, of the incident beam particles. It uses two auxiliary files containing
    masses and stopping powers (/usr/local/tab/massa.tab /usr/local/tab/scoef.88).

    The program asks the reaction, to be specified by 1) beam particle, 2)
    target nuclide, 3) detected particle, and 4) residual nucleus.
    If the residual nucleus is not in its ground state, its excitation energy
    may be given (level of particles 4).

    To obtain the depth energy relation of the beam particles, give a 0 instead
    of particle 2; particles 3 and 4 are then not asked.

    The target and detector orientation are specified by angles phi_in, phi_out
    and theta. The 'depth' in the output is incoming path length * cos(phi_in).

    Next, the stopping material is to be given. One may specify either a single
    element (e.g. Si), or a compound (e.g. Si1O2). Only in the second case,
    i.e. for compunds, the density is also asked.

    Then the beam energy (MeV) is asked, and the channel number corresponding
    to the maximum energy of particles 3. Instead of this channel number one
    may give the calibration in keV per channel. (Input is then a negative
    value, -keV/channel.) The energy of channel 0, is also asked.

    Next, select one of the output options, D, E, or N. For option D, the depth
    (in nm) corresponding to each channel is stored. Additional input is the
    name of the output array. This should be a 'double' or 'float' array of
    size at least 1024. In the second case (option E) the depth-energy relation
    is stored in an array, assumed to have 1024 channels. The channel numbers
    of this array correspond to depth, and the channel contents to energy in keV.
    Channel 0 corresponds to depth 0. In this case the program also asks for the
    number of nm per channel.

    Further the program asks for an energy step (MeV) to be used as interval
    in the energy of particles 1, in the calculations. A default value is
    obtained by just giving a comma at this point.

    After the calculations are finished, the program outputs the number of
    steps actually used, and one may decide to repeat the calculation with a
    smaller energy step if this number is too low. (the maximum nr of steps the
    program can handle is 128, if the number used is less then about 30 results
    may be unreliable).

    After establishing the depth-energy relation the program again asks for
    input, in the form of a list of numbers, interpreted as follows: positive
    numbers are depth in nm. negative numbers between 0 and -10 are (minus)
    energies in MeV. negative numbers less then -10 are (minus) channels
    numbers.

    For each of these the depth and corresponding energy are output. To quit
    command DEPth, give a ';'

ISOtopes

    ISO <list of element names>

      Scans the file 'massa.tab' and list all isotopes of the given elements.
      The quantities Z,N,A,and mass excess are output.
      In addition the average mass and density for the natural composition are
      listed, as found from file 'scoef.88'.
      Leave the command ISO by giving a `;' or an empty line.

      Example: ISO LI,F;

MASS

    MASS <projectile> <energy-in> <scattering angle> <energy-out>

      Calculates the mass of the scattering nucleus, assuming elastic
      scattering. Relativistically correct.

      projectile : mnemonic for the projectile (e.g. P,A,22NE)
      energy-in  : kinetic energy in MeV
      scattering angle : in degrees (forward is 0)
      energy-out : kinetic energy in MeV after the collision

      The mass is printed in atomic mass units.

Spectrum simulation

RBSIM

    This commands builds a simulated energy spectrum of particles from a
    nuclear reaction, taking into account:
    1. Y(z) : the concentration and/or nuclear encounter probability as a
              function of depth. Read from an ASCII file or an -.FLX file.
    2. The energy loss and straggling in the incoming trajectory, using
       Ziegler et al.'s stopping powers and Bohr straggling, or, read in
       from an -.FLX file.
    3. The energy loss and straggling in the outgoing trajectory, using
       Ziegler et al.'s stopping powers and Bohr straggling.
    4. X(E): the reaction cross section as a function of energy.
       Options:
         0 : uniform, or already included in Y(z).
         1 : Rutherford cross section
         2 : read in from an ASCII file
    5. The detector resolution (FWHM).
    Restriction. An approximation is made that is only valid if the yield
    varies smoothly with energy, i.e. at any depth the cross section is
    assumed constant, although the energy has a distribution of finite width.

    Input:
      a. FLX file name (give a comma or <return> if not needed)
      b. ASCII file name (give a comma or <return> if not needed)
      c. Random (=uniform) Y(z)? (type Yes or No)
      d. Random energy loss? (Yes or No;
                              if 'No': use simulated energy loss of FLX file)
      e. Random straggling? (ditto)
      f. Cross section type (0,1,2; see above)
      g. The following input depends on the options chosen.
           If a cross section of type 1 was specified the program asks the
           subfile name in the ASCII file.
           If a non-random Y(z) was specified, the program asks what to read:
             1=ASCII file
             2=FLX file (FLUX V2.0)
             3=FLX file (FLUX V3.0)
           After option 1, a subfile name is asked, and also the value of the
           beam energy.
           After option 2, a sequence number is asked.
           After option 3, a sequence number is asked, and also the atom
           species (1 or 2) that is involved in the reaction.
      If a random Y(z) was specified, the program asks for the beam energy
      and also of the maximum depth in nm.
      h. Specify the reaction, example: 4HE 28SI 4HE 28SI
      i. phi-in, phi-out, theta
      j. Substrate composition (example SI1O2)
      k. Density (g/cm3), average atomic number A
         If the substrate consists of a single element, the density and A are
         nor asked, but looked up instead.
      l. Thickness of extra surface layer (nm)
         Fraction of bulk atoms in surface layer.
         (this is to "add" an extra surface peak to the Y(z) specified
         it does not take the different stopping into account)
      m. Spectrum calibration (+ch where E=Emax, or, -keV/channel)
         zero crossing (E(keV) corresponding to ch 0)
      n. Spectrum name (spectrum should be 1024 channels or more)
         optional spectrum for depth scale
      o. Detector resolution (keV FWHM)
      p. Fin,Fuit: multiplicative factors for the stopping powers

Fit, smooth, interpolate

LLSQ

    input parameters                   comment
    ----------------                   -------
    llsq                               Linear Least SQuares
    array of x values                  all arrays must be of type float
    array to be fitted                 data
    array of *squared* errors          may be the same as the data array
    array for fitted values
    first channel number to be fitted  the same for all 4 arrays
    last channel number to be fitted
    Option `L' or `P'                  P or p for fitting power series
                                       L or l for Legendre polynomials.
    NR. of L-values.                   number of fit parameters
    L-values                           example: for quadratic curve, type
                                       option = P,
                                       NR. of L-values = 3
                                       L-values = 0, 1, 2
    Option `Y' or `N' or `S'           After the fit, type
                                       S or s to quit
                                       Y or N to repeat with different
                                       NR. of L-values or different
                                       L-values
                                       
    The output contains the coefficients A(L) in the specified series:
      L=Legendre polynomials : Y=sum{A(L) × PL(cos(x)}
      P=Power series :         Y=sum{A(L) × x**L}

SPLINE

    spline <option> <further parameters>

      <option> = one of the characters s,0,1,2,3
      option s : smooth a given set of data, using fit with cubic splines.
      option n=0,1,2,3 : calculate n'th derivative from a previous spline fit.

      input for option s :
        A : array to be smoothed
        B : array for result; default is array A.
        C : array of relative or absolute squared errors; default is array A.
        f, l : first and last channel to be fitted; default is the whole array A
        variance parameter : 1 if absolute squared errors were given,
                             see also further on.
        A number of CIO arrays is created with names splinx, spliny, ...
        The arrays splinx, spliny, and splinc are used by options 0-3,
        and should not be modified between spline s and spline 0,1,2, or 3.

      input for options 0,1,2,3 :
        xarray : an array of x values; default is splinx
        output : an array of results.
        f, l :   first and last channel containing x values; default is the
                 whole xarray

        For every x value the interpolated function value (option 0) or the
        value of the n'th derivative (options 1,2,3) is stored in the
        corresponding channel of the result array.

C   AUTHOR              - M.F.HUTCHINSON
C                         CSIRO DIVISION OF MATHEMATICS AND STATISTICS
C                         P.O. BOX 1965
C                         CANBERRA, ACT 2601
C                         AUSTRALIA
C
C   LATEST REVISION     - 15 AUGUST 1985
C
C   PURPOSE             - CUBIC SPLINE DATA SMOOTHER
C

 1.  Craven, Peter and Wahba, Grace, "Smoothing noisy data with spline
     functions", Numer. Math. 31, 377-403 (1979).
 2.  Hutchinson, M.F. and de Hoog, F.R., "Smoothing noisy data with spline
     functions", Numer. Math. (in press).
 3.  Reinsch, C.H., "Smoothing by spline functions", Numer. Math. 10,
     177-183 (1967).
 4.  Wahba, Grace, "Bayesian 'confidence intervals' for the cross-validated
     smoothing spline", J.R.Statist. Soc. B 45, 133-150 (1983).

EXAMPLE

A sequence of 50 data points are generated by adding a random variable with
uniform density in the interval [-0.3,0.3] to the curve y=sin(3*pi*x/2).
The abscissae are unequally spaced in [0,1]. Point standard error estimates
are returned in SE by setting JOB to 1. The error variance estimate is
returned in VAR. It compares favourably with the true value of 0.03. The
IMSL function GGUBFS is used to generate sample values of a uniform
variable on [0,1].

INPUT:

      INTEGER          N,IC,JOB,IER
      DOUBLE PRECISION X(50),F(50),Y(50),DF(50),C(49,3),WK(400),
     *                 VAR,SE(50)
      DOUBLE PRECISION GGUBFS,DSEED,DN
      DATA DSEED/1.2345D4/
 C
      N=50
      IC=49
      JOB=1
      VAR=-1.0D0
      DN=N
 C
      DO 10 I=1,N
      X(I)=(I - 0.5)/DN + (2.0*GGUBFS(DSEED) - 1.0)/(3.0*DN)
      F(I)=DSIN(4.71238*X(I)) + (2.0*GGUBFS(DSEED) - 1.0)*0.3
      DF(I)=1.0D0
  10  CONTINUE
      CALL CUBGCV(X,F,DF,N,Y,C,IC,VAR,JOB,SE,WK,IER)
       .
       .
       .
      END

OUTPUT:

 IER = 0
 VAR = 0.0279
 GENERALIZED CROSS VALIDATION = 0.0318
 MEAN SQUARE RESIDUAL = 0.0246
 RESIDUAL DEGREES OF FREEDOM = 43.97
 FOR CHECKING PURPOSES THE FOLLOWING OUTPUT IS GIVEN:

 X(1)  = 0.0046    F(1)  =  0.2222     Y(1)  =  0.0342     SE(1)  = 0.1004
 X(21) = 0.4126    F(21) =  0.9810     Y(21) =  0.9069     SE(21) = 0.0525
 X(41) = 0.8087    F(41) = -0.5392     Y(41) = -0.6242     SE(41) = 0.0541

SPLINT

    Interpolation of a function using cubic splines.

    a) setup for spline interpolation :

       splint s <degree> <x-array> <y-array> <first> <last> [derivs] 0

         <degree>            : 3 or 5, degree of the spline.
         <x-array> <y-array> : arrays containing points where the function
                               to be interpolated has prescribed values
         <x-array>           : the x-values in increasing order
         <y-array>           : the function values.
         <first> <last>      :     indices of the range in <x-array> <y-array>
         [derivs]            : It is possible to prescribe the value of derivatives
                               at any of the points given above. To do so enter one
                               or more lines, containing
                                 <n> <ch> <value>
                                 where <n> = 1,2 or 3 for 1'st, 2'nd, 3'rd derivative
                                       <ch> = index in <x-array> <y-array>
                                       <value> = value of the derivative
         0                   : the input is terminated by <n> = 0, or a semicolon

    b) interpolation

       splint <id> <x-array> <y-array> <first> <last>

         <id>           : has possible values 0, 1,2,3 and specifies whether the
                          function itself, or one of its derivatives is wanted
         <x-array>      : contains x values where the function should be caluclated.
         <y-array>      : array where the result is return in
         <first> <last> : indicate the range in <x-array> <y-array>

 This command uses a Fortran procedure from netlib:
 C     ALGORITHM 485 COLLECTED ALGORITHMS FROM ACM.
 C     ALGORITHM APPEARED IN COMM. ACM, VOL. 17, NO. 09,
 C     P. 526.
       SUBROUTINE GSF(N, M, X, Y, Z, IM, C, IDET)
       INTEGER N, M
       REAL X(N), Y(4,N)
       INTEGER Z(N)
       INTEGER IM(4,N)
       REAL C(8,N)
       INTEGER IDET

 C INPUT   N,M,X,Z,IM,Y--
 C N IS A POSITIVE INTEGER GIVING THE NUMBER OF KNOTS
 C M IS A POSITIVE INTEGER DETERMINING THE DEGREE
 C 2*M-1 OF THE SPLINE
 C X IS AN ARRAY OF REAL NUMBERS WITH
 C     X(1).LT.X(2).LT.....LT.X(N)
 C Z IS AN ARRAY OF INTEGERS SUCH THAT
 C     0.LT.Z(I).LE.M, I=1,2,...,N
 C IM IS AN INTEGER ARRAY WITH
 C     1.LE.M(1,J).LT.....LT.IM(Z(J),J).LE.M,
 C J=1,2,...,N.
 C Y IS AN ARRAY OF REAL NUMBERS
 C OUTPUT  C,IDET--
 C THE COLUMN VECTOR C(1,J),...,C(2*M,J) CONTAINS THE
 C COEFFICIENTS OF THE SPLINE IN THE INTERVAL X(J-1)
 C TO X(J). IDET IS SET TO ZERO IF A SINGULAR SYSTEM
 C IS ENCOUNTERED OTHERWISE, IDET IS 1.
 C THE SUBROUTINE  GSF(N,M,X,Y,Z,IM,C,IDET)  COMPUTES
 C THE COEFFICIENTS OF THE INTERPOLATING G-SPLINE. THE
 C PARAMETERS N,M,X,Y,Z, AND  IM  ARE INPUT. N  AND  M
 C ARE THE POSITIVE INTEGERS OF SECTION 1  WHICH GIVE
 C THE NUMBER OF X S AND DETERMINE THE DEGREE OF THE
 C SPLINE, RESPECTIVELY. X  MUST BE AN ARRAY OF REAL
 C NUMBERS WITH X(1).LT.X(2).LT.....LT.X(N)  AND Z  IS
 C AN ARRAY OF POSITIVE INTEGERS NONE OF WHICH SHOULD
 C EXCEED M. X CONTAINS THE POINTS WHERE HB-DATA IS
 C PRESCRIBED AND  Z DESCRIBES THE NUMBER OF PIECES OF
 C DATA AT EACH SUCH POINT. IM  IS AN INTEGER ARRAY
 C WITH
 C  1.LE.IM(1,J).LT.IM(2,J).LT.....LT.IM(Z(J),L).LE.M,
 C              J=1,2,...,N.
 C THE J TH COLUMN OF  IM  IS A LIST OF WHICH
 C DERIVATIVES (SHIFTED UP BY 1) ARE SPECIFIED AT
 C X(J).  THE DATA FOR THE HB-INTERPOLATION PROBLEM IS
 C ENTERED IN THE ARRAY  Y. Y(I,J) SHOULD CONTAIN THE
 C VALUE ASSIGNED TO THE IM(I,J)-1 ST DERIVATIVE OF
 C THE SPLINE EVALUATED AT  X(J). THE PARAMETERS  C
 C AND  IDET  ARE OUTPUT OF GSF. AFTER EXECUTION, THE
 C ARRAY C WILL CONTAIN THE COEFFICIENTS OF THE
 C SPLINE. IN PARTICULAR, THE COLUMN VECTOR
 C C(1,J),...,C(2*M,J) CONTAINS THE COEFFICIENTS OF
 C THE POLYNOMIALS P(J)  DESCRIBED BY EQUATIONS (6),
 C J=1,2,...,N+1. SUBROUTINE GSF CALLS ON SUBROUTINE
 C INTCON,SMOCON, AND TRISYS WHICH MUST BE LOADED WITH
 C THE MAIN PROGRAM.

 [ ... ]

       FUNCTION GVAL(T, ID, N, M, X, C)
       REAL T
       INTEGER ID
       INTEGER N, M
       REAL X(N)
       REAL C(8,N)

 C INPUT T,ID,N,M,X,C
 C THE PARAMETERS N,M,X,C ARE AS IN GSF  AND
 C COMPLETELY DESCRIBE THE G-SPLINE.
 C T IS A REAL NUMBER AND ID A POSITIVE INTEGER.
 C GVAL PRODUCES THE ID-1 ST DERIVATIVE OF THE SPLINE
 C AT T. GVAL AUTOMATICALLY PRODUCES 0 IF ID.GT.M*2