#! /bin/sh
# a script to set up the interface code for accessing
# FormCalc-generated code in Mathematica
# this file is part of FormCalc
# last modified 14 Jul 10 th


if [ $# -ne 1 ] ; then
  echo "Usage:  $0 file.F   -- generates file.tm for file.F" 1>&2
  exit 1
fi

file="$1"
base=`basename "$file" .F`
name=`echo "$base" | sed s/[-_]//`

pol=`awk '/^#include *".to.\.F"/ {
  print substr("UUUUUUUUUU", 1, int(substr($2,2,1)) + int(substr($2,5,1)));
  exit(0);
}' process.h`

set -- `sed '
s/.*call MmaGetInteger(/i /
s/.*call MmaGetIntegerList(/il /
s/.*call MmaGetReal(/r /
s/.*call MmaGetRealList(/rl /
s/.*call MmaGetComplex(/c /
s/.*call MmaGetComplexList(/cl /
y/_,)/$  /
t
d' $file`

while [ $# -gt 0 ] ; do
  case $1 in
  i)	argsQ="$argsQ, $2_Integer"
	argsL="$argsL, $2_"
	argsR="$argsR, $2"
	argsN="$argsN, $2"
	argsT="$argsT, Integer"
	argsC="$argsC, cint $2"
	argsS="$argsS, $2"
	;;
  il)	argsQ="$argsQ, $2_?(il[$3])"
	argsL="$argsL, $2_"
	argsR="$argsR, $2"
	argsN="$argsN, $2"
	argsT="$argsT, IntegerList"
	argsC="$argsC, int *$2, long n_$2"
	argsS="$argsS, $2, n_$2"
	shift
	;;
  r)	argsQ="$argsQ, $2_?r"
	argsL="$argsL, $2_"
	argsR="$argsR, $2"
	argsN="$argsN, N[$2]"
	argsT="$argsT, Real"
	argsC="$argsC, creal $2"
	argsS="$argsS, $2"
	;;
  rl)	argsQ="$argsQ, $2_?(rl[$3])"
	argsL="$argsL, $2_"
	argsR="$argsR, $2"
	argsN="$argsN, N[$2]"
	argsT="$argsT, RealList"
	argsC="$argsC, real *$2, long n_$2"
	argsS="$argsS, $2, n_$2"
	shift
	;;
  c)	argsQ="$argsQ, $2_?c"
	argsL="$argsL, $2_"
	argsR="$argsR, $2"
	argsN="$argsN, N[Re[$2]], N[Im[$2]]"
	argsT="$argsT, Real, Real"
	argsC="$argsC, creal re$2, creal im$2"
	argsS="$argsS, (cplx){re$2, im$2}"
	;;
  cl)	argsQ="$argsQ, $2_?(cl[$3])"
	argsL="$argsL, $2_"
	argsR="$argsR, $2"
	argsN="$argsN, N[Flatten[{Re[#], Im[#]}&/@ $2]]"
	argsT="$argsT, RealList"
	argsC="$argsC, real *$2, long n_$2"
	argsS="$argsS, $2, n_$2"
	shift
	;;
  esac
  shift 2
done

file="`dirname $file`/$base.tm"

cat << _EOF_ > $file
:Evaluate: BeginPackage["XSection\`"]

:Evaluate: Options[$name] = {
  Polarizations -> "$pol",
  Serial -> {},
  SetNumber -> 1,
  ParaHead -> Para,
  DataHead -> Data,
  LogFile -> "" }

:Evaluate: DataRow::usage = "DataRow[v, r, e] stores a row of data from FormCalc-generated Fortran code.
	The values of the independent variables are given in v, the integration results in r, and the corresponding integration errors in e."

:Evaluate: Begin["\`$name\`"]

:Evaluate: r = Head[# + 1.] === Real &

:Evaluate: c = Head[# + 1. I] === Complex &

:Evaluate: il[n_] = Length[#] === n && VectorQ[#, IntegerQ] &

:Evaluate: rl[n_] = Length[#] === n && VectorQ[#, r] &

:Evaluate: cl[n_] = Length[#] === n && VectorQ[#, c] &

:Evaluate: tocl = Chop[Apply[Complex, Partition[#, 2], 1], 2^-51] &

:Evaluate: $name[sqrtS_$argsQ, opt___Rule] :=
  Block[ {pol, serial, setno, phead, dhead, file},
    {pol, serial, setno, phead, dhead, file} =
      {Polarizations, Serial, SetNumber, ParaHead, DataHead, LogFile} /.
        {opt} /. Options[$name];
    m$name[setno, ToString[phead], ToString[dhead],
      ToString[file], ToString[pol],
      Flatten[{sqrtS}], Flatten[{serial}]$argsR]
  ]

:Begin:
:Function: m$name
:Pattern: m$name[setno_Integer, phead_String, dhead_String,
  file_String, pol_String,
  {sqrtSfrom_, sqrtSto_:-1, sqrtSstep_:10},
  {serialfrom_:0, serialto_:2^30, serialstep_:1}$argsL]
:Arguments: {setno, phead, dhead,
  file, pol,
  N[sqrtSfrom], N[sqrtSto], N[sqrtSstep],
  serialfrom, serialto, serialstep$argsN}
:ArgumentTypes: {Integer, String, String,
  String, String,
  Real, Real, Real,
  Integer, Integer, Integer$argsT}
:ReturnType: Manual
:End:

:Evaluate: _m$name := (Message[$name::args]; Abort[])

:Evaluate: $name::args = "Invalid arguments used in $name."

:Evaluate: $name::redir = "Cannot redirect output to \`\`."

:Evaluate: End[]

:Evaluate: EndPackage[]


/*
	$base.tm
		Mathematica interface to FormCalc-generated code
		generated by `basename $0`
		`date`
*/

#include "mathlink.h"

#ifndef MLCONST
#define MLCONST const
#endif

#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <fcntl.h>

typedef MLCONST char *string;
typedef const int cint;
typedef double real;
typedef const real creal;
typedef struct { double re, im; } cplx;
typedef const cplx ccplx;

#if UNDERSCORE
#define FORTRAN(s) s##_
#else
#define FORTRAN(s) s
#endif

extern void FORTRAN(processini)(int *fail, string pol,
  creal *sqrtSfrom, creal *sqrtSto, creal *sqrtSstep, cint pollen);

extern void FORTRAN(parameterscan)(string dir,
  cint *serialfrom, cint *serialto, cint *serialstep, cint dirlen);

extern void FORTRAN(fortranflush)(void);

static va_list xsection_ap;
static string parahead, datahead;
static int setnumber;

/******************************************************************/

static inline void MLPutComplex(MLINK mlp, ccplx c)
{
  if( c.im == 0 ) MLPutReal(mlp, c.re);
  else {
    MLPutFunction(mlp, "Complex", 2);
    MLPutReal(mlp, c.re);
    MLPutReal(mlp, c.im);
  }
}

/******************************************************************/

static inline void MLSendPacket(MLINK mlp)
{
  int pkt;

  MLEndPacket(mlp);

  do {
    pkt = MLNextPacket(mlp);
    MLNewPacket(mlp);
  } while( pkt != RETURNPKT );
}

/******************************************************************/

static inline void MLEmitMessage(MLINK mlp, string tag, string arg)
{
  MLPutFunction(mlp, "EvaluatePacket", 1);

  MLPutFunction(mlp, "Message", 2);
  MLPutFunction(mlp, "MessageName", 2);
  MLPutSymbol(mlp, "$name");
  MLPutString(mlp, tag);
  MLPutString(mlp, arg);

  MLSendPacket(mlp);
}

/******************************************************************/

static void MLPutRule(MLINK mlp, string s, cint len)
{
  MLPutFunction(mlp, "Rule", 2);
  MLPutFunction(mlp, "ToExpression", 1);
  MLPutByteString(mlp, (unsigned MLCONST char *)s, len);

/*
  char *const tmp = malloc(len + 1);
  char *d;
  int i;

  if( tmp == NULL ) exit(1);
  d = tmp;

  for( i = 0; i < len; ++i ) {
    char c = *s++;
    switch( c ) {
    case '_':
      c = '$';
      break;
    case '(':
      c = '[';
      break;
    case ')':
      c = ']';
      break;
    }
    *d++ = c;
  }
  *d++ = 0;

  MLPutFunction(mlp, "Rule", 2);
  MLPutFunction(mlp, "ToExpression", 1);
  MLPutString(mlp, tmp);

  free(tmp);
*/
}

/******************************************************************/

static void v$name(
  cint serialfrom, cint serialto, cint serialstep, ...)
{
  va_start(xsection_ap, serialstep);
  FORTRAN(parameterscan)("", &serialfrom, &serialto, &serialstep, 0);
  va_end(xsection_ap);
}

/******************************************************************/

static void m$name(cint setno, string phead, string dhead,
  string file, string pol,
  creal sqrtSfrom, creal sqrtSto, creal sqrtSstep,
  cint serialfrom, cint serialto, cint serialstep$argsC)
{
  int fail;
  cint prevstdout = dup(1);

  if( *file == 0 ) dup2(2, 1);
  else {
    cint logfile = creat(file, 0644);
    if( logfile == -1 ) {
      MLEmitMessage(stdlink, "redir", file);
      MLPutSymbol(stdlink, "$Failed");
      MLEndPacket(stdlink);
      return;
    }
    dup2(logfile, 1);
    close(logfile);
  }

  setnumber = setno - 1;
  parahead = phead;
  datahead = dhead;

  FORTRAN(processini)(&fail, pol,
    &sqrtSfrom, &sqrtSto, &sqrtSstep, strlen(pol));
  v$name(serialfrom, serialto, serialstep$argsS);

  FORTRAN(fortranflush)();
  dup2(prevstdout, 1);
  close(prevstdout);

  MLPutInteger(stdlink, setnumber);
  MLEndPacket(stdlink);
}

/******************************************************************/

void FORTRAN(mmagetinteger)(int *i)
{
  *i = va_arg(xsection_ap, int);
}

/******************************************************************/

void FORTRAN(mmagetintegerlist)(int *i, cint *n)
{
  const int *l = va_arg(xsection_ap, int *);
  const long nl = va_arg(xsection_ap, long);
  memcpy(i, l, nl*sizeof(int));
	/* note: no buffer overflow due to pattern in Mma */
}

/******************************************************************/

void FORTRAN(mmagetreal)(real *r)
{
  *r = va_arg(xsection_ap, real);
}

/******************************************************************/

void FORTRAN(mmagetreallist)(real *r, cint *n)
{
  const real *l = va_arg(xsection_ap, real *);
  const long nl = va_arg(xsection_ap, long);
  memcpy(r, l, nl*sizeof(real));
}

/******************************************************************/

void FORTRAN(mmagetcomplex)(cplx *c)
{
  c->re = va_arg(xsection_ap, real);
  c->im = va_arg(xsection_ap, real);
}

/******************************************************************/

void FORTRAN(mmagetcomplexlist)(cplx *c, cint *n)
{
  const real *l = va_arg(xsection_ap, real *);
  const long nl = va_arg(xsection_ap, long);
  memcpy(c, l, nl*sizeof(real));
}

/******************************************************************/

void FORTRAN(mmasetpara)()
{
  MLPutFunction(stdlink, "EvaluatePacket", 1);

  MLPutFunction(stdlink, "Set", 2);
  MLPutFunction(stdlink, parahead, 1);
  MLPutInteger(stdlink, ++setnumber);

  MLPutFunction(stdlink, "List", 1);
}

/******************************************************************/

void FORTRAN(mmasetdata)()
{
  MLPutFunction(stdlink, "EvaluatePacket", 1);

  MLPutFunction(stdlink, "Set", 2);
  MLPutFunction(stdlink, datahead, 1);
  MLPutInteger(stdlink, setnumber);

  MLPutFunction(stdlink, "List", 1);
}

/******************************************************************/

void FORTRAN(mmaendset)()
{
  MLPutFunction(stdlink, "Sequence", 0);
  MLSendPacket(stdlink);
}

/******************************************************************/

void FORTRAN(mmaputinteger)(string name, cint *i, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutRule(stdlink, name, len);
  MLPutInteger(stdlink, *i);
}

/******************************************************************/

void FORTRAN(mmaputintegerlist)(string name, int *i, cint *n, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutRule(stdlink, name, len);
  MLPutIntegerList(stdlink, i, *n);
}

/******************************************************************/

void FORTRAN(mmaputreal)(string name, creal *r, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutRule(stdlink, name, len);
  MLPutReal(stdlink, *r);
}

/******************************************************************/

void FORTRAN(mmaputreallist)(string name, real *r, cint *n, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutRule(stdlink, name, len);
  MLPutRealList(stdlink, r, *n);
}

/******************************************************************/

void FORTRAN(mmaputcomplex)(string name, ccplx *c, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutRule(stdlink, name, len);
  MLPutComplex(stdlink, *c);
}

/******************************************************************/

void FORTRAN(mmaputcomplexlist)(string name, cplx *c, cint *n, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutRule(stdlink, name, len);
  MLPutFunction(stdlink, "XSection\`$name\`tocl", 1);
  MLPutRealList(stdlink, (real *)c, *n*2);
}

/******************************************************************/

void FORTRAN(mmadata)(real *show, cint *nshow,
  real *result, real *error, cint *ncomp)
{
  MLPutFunction(stdlink, "Sequence", 2);

  MLPutFunction(stdlink, "DataRow", 3);
  MLPutRealList(stdlink, show, *nshow);
  MLPutRealList(stdlink, result, *ncomp);
  MLPutRealList(stdlink, error, *ncomp);
}

/******************************************************************/

int main(int argc, char **argv)
{
  return MLMain(argc, argv);
}

_EOF_

