#! /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 23 Apr 13 th


test $# -ne 1 && {
  echo "Usage:  FC FFLAGS file.F | $0 file.tm   -- generates file.tm for file.F" 1>&2
  exit 1
}

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

prep=`cat -`

pol=`echo "$prep" | awk '/len_trim\(pol\)/ { print substr("UUUUUUUUUU", 1, $4); }'`

set -- `echo "$prep" | 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'`

while test $# -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, cRealType $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, RealType *$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, cRealType re$2, cRealType im$2"
	argsS="$argsS, (ComplexType){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, RealType *$2, long n_$2"
	argsS="$argsS, $2, n_$2"
	shift
	;;
  esac
  shift 2
done

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

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

: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, args, serial, setno, phead, dhead, logfile, pset, dset},
    {pol, args, serial, setno, phead, dhead, logfile} =
      {Polarizations, Args, Serial, SetNumber,
        ParaHead, DataHead, LogFile} /. {opt} /. Options[$name];
    args = StringJoin[{ToString[#], "\000"}&/@ args];
    pset = With[{h = phead}, (h[#1] = #2)&];
    dset = With[{h = dhead}, (h[#1] = #2)&];
    m$name[setno, ToString[logfile], ToString[pol], args,
      Flatten[{sqrtS}], Flatten[{serial}]$argsR]
  ]

:Begin:
:Function: m$name
:Pattern: m$name[setno_Integer, logfile_String,
  pol_String, args_String,
  {sqrtSfrom_, sqrtSto_:-1, sqrtSstep_:10},
  {serialfrom_:0, serialto_:2^30, serialstep_:1}$argsL]
:Arguments: {setno, logfile, pol, args,
  N[sqrtSfrom], N[sqrtSto], N[sqrtSstep],
  serialfrom, serialto, serialstep$argsN}
:ArgumentTypes: {Integer, ByteString,
  ByteString, ByteString,
  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>
#include <pthread.h>

typedef unsigned char byte;
typedef MLCONST byte *string;
typedef byte argstr[128];
typedef const int cint;
typedef double RealType;
typedef const RealType cRealType;
typedef struct { double re, im; } ComplexType;
typedef const ComplexType cComplexType;

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

#ifdef __cplusplus
extern "C" {
#endif

extern void FORTRAN(processini)(int *fail, string pol,
  cRealType *sqrtSfrom, cRealType *sqrtSto, cRealType *sqrtSstep,
  cint *argc, argstr argv[],
  cint pollen, cint argstrlen);

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

extern void FORTRAN(processexi)(void);

extern void FORTRAN(fortranflush)(void);

void FORTRAN(mmagetinteger)(int *i);
void FORTRAN(mmagetintegerlist)(int *i, cint *n);
void FORTRAN(mmagetreal)(RealType *r);
void FORTRAN(mmagetreallist)(RealType *r, cint *n);
void FORTRAN(mmagetcomplex)(ComplexType *c);
void FORTRAN(mmagetcomplexlist)(ComplexType *c, cint *n);

void FORTRAN(mmaputinteger)(string name, cint *i, cint len);
void FORTRAN(mmaputintegerlist)(string name, int *i, cint *n, cint len);
void FORTRAN(mmaputreal)(string name, cRealType *r, cint len);
void FORTRAN(mmaputreallist)(string name, RealType *r, cint *n, cint len);
void FORTRAN(mmaputcomplex)(string name, cComplexType *c, cint len);
void FORTRAN(mmaputcomplexlist)(string name, ComplexType *c, cint *n, cint len);

void FORTRAN(mmasetpara)(void);
void FORTRAN(mmasetdata)(void);
void FORTRAN(mmaendset)(void);
void FORTRAN(mmadata)(RealType *show, cint *nshow,
  RealType *result, RealType *error, cint *ncomp);

#ifdef __cplusplus
}
#endif

static va_list xsection_ap;
static int setnumber;
static int stdoutorig;
static pthread_mutex_t mmamutex;

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

#define Strlen(s) strlen((const char *)s)
#define Strncpy(s1, s2, n) strncpy((char *)s1, (const char *)s2, n)

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

static inline void MLSendPacket(MLINK mlp)
{
  MLEndPacket(mlp);
  MLNextPacket(mlp);
  MLNewPacket(mlp);
}

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

static inline void MLPrint(MLINK mlp, string s, cint len)
{
  pthread_mutex_lock(&mmamutex);
  MLPutFunction(mlp, "EvaluatePacket", 1);
  MLPutFunction(mlp, "WriteString", 2);
  MLPutString(mlp, "stdout");
  MLPutByteString(mlp, s, len);
  MLSendPacket(mlp);
  pthread_mutex_unlock(&mmamutex);
}

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

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

/*
  byte tmp[len + 1], *d = tmp;
  int i;

  for( i = 0; i < len; ++i ) {
    byte 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);
  MLPutByteString(mlp, tmp, d - tmp);

  free(tmp);
*/
}

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

static void *MLstdout(void *fd)
{
  byte s[1024];
  ssize_t n;

  while( (n = read(*(int *)fd, s, sizeof s)) > 0 )
    MLPrint(stdlink, s, n);
  return NULL;
}

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

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

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

static void m$name(cint setno, string logfile, cint logfilelen,
  string pol, cint pollen, string args, cint argslen,
  cRealType sqrtSfrom, cRealType sqrtSto, cRealType sqrtSstep,
  cint serialfrom, cint serialto, cint serialstep$argsC)
{
  int fail, stdoutpipe[2], argc;
  enum { argmax = 10 };
  argstr argv[argmax];
  string argend;
  pthread_t tid = 0;
  void *dummy;
  char file[logfilelen + 1];

  memcpy(file, logfile, logfilelen);
  file[logfilelen] = 0;
  if( strcmp(file, "stdout") == 0 ) {
    if( pipe(stdoutpipe) == -1 ||
        pthread_create(&tid, NULL, MLstdout, stdoutpipe) != 0 )
      stdoutpipe[1] = -1;
  }
  else stdoutpipe[1] = creat(file, 0644);

  if( stdoutpipe[1] == -1 ) {
    MLPutFunction(stdlink, "EvaluatePacket", 1);
    MLPutFunction(stdlink, "Message", 2);
    MLPutFunction(stdlink, "MessageName", 2);
    MLPutSymbol(stdlink, "$name");
    MLPutString(stdlink, "redir");
    MLPutByteString(stdlink, logfile, logfilelen);
    MLSendPacket(stdlink);

    MLPutSymbol(stdlink, "$Failed");
    MLEndPacket(stdlink);
    return;
  }

  pthread_mutex_init(&mmamutex, NULL);

  dup2(stdoutpipe[1], 1);
  close(stdoutpipe[1]);

  setnumber = setno - 1;

  argend = args + argslen;
  argc = 0;
  while( argc < argmax && args < argend ) {
    Strncpy(argv[argc++], args, sizeof(argstr));
    args = (byte *)memchr(args, 0, argend - args) + 1;
  }

  FORTRAN(processini)(&fail, pol,
    &sqrtSfrom, &sqrtSto, &sqrtSstep,
    &argc, argv,
    pollen, sizeof(argstr));
  v$name(serialfrom, serialto, serialstep$argsS);
  FORTRAN(processexi)();

  FORTRAN(fortranflush)();
  fflush(stdout);
  dup2(stdoutorig, 1);
  if( tid ) pthread_join(tid, &dummy);
  pthread_mutex_destroy(&mmamutex);

  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)(RealType *r)
{
  *r = va_arg(xsection_ap, RealType);
}

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

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

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

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

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

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

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

void FORTRAN(mmasetpara)()
{
  pthread_mutex_lock(&mmamutex);
  MLPutFunction(stdlink, "EvaluatePacket", 1);
  MLPutFunction(stdlink, "XSection\`$name\`pset", 2);
  MLPutInteger(stdlink, ++setnumber);
  MLPutFunction(stdlink, "List", 1);
}

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

void FORTRAN(mmasetdata)()
{
  pthread_mutex_lock(&mmamutex);
  MLPutFunction(stdlink, "EvaluatePacket", 1);
  MLPutFunction(stdlink, "XSection\`$name\`dset", 2);
  MLPutInteger(stdlink, setnumber);
  MLPutFunction(stdlink, "List", 1);
}

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

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

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

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, cRealType *r, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);
  MLPutRule(stdlink, name, len);
  MLPutReal(stdlink, *r);
}

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

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

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

void FORTRAN(mmaputcomplex)(string name, cComplexType *c, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);
  MLPutRule(stdlink, name, len);
  if( c->im == 0 ) MLPutReal(stdlink, c->re);
  else {
    MLPutFunction(stdlink, "Complex", 2);
    MLPutReal(stdlink, c->re);
    MLPutReal(stdlink, c->im);
  }
}

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

void FORTRAN(mmaputcomplexlist)(string name, ComplexType *c, cint *n, cint len)
{
  MLPutFunction(stdlink, "Sequence", 2);
  MLPutRule(stdlink, name, len);
  MLPutFunction(stdlink, "XSection\`$name\`tocl", 1);
  MLPutRealList(stdlink, (RealType *)c, *n*2);
}

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

void FORTRAN(mmadata)(RealType *show, cint *nshow,
  RealType *result, RealType *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)
{
  int fd;

	/* make sure a pipe will not overlap with 0, 1, 2 */
  do fd = open("/dev/null", O_WRONLY); while( fd <= 2 );
  close(fd);
  stdoutorig = dup(1);

  return MLMain(argc, argv);
}

_EOF_

