(*^

::[paletteColors = 128; 
	fontset = title, "New York", 24, L3, center, bold, nohscroll;
	fontset = subtitle, "New York", 18, L2, center, bold, nohscroll;
	fontset = subsubtitle, "New York", 14, L2, center, bold, nohscroll;
	fontset = section, "New York", 14, L2, bold, nohscroll, grayBox;
	fontset = subsection, "New York", 12, L2, bold, nohscroll, blackBox;
	fontset = subsubsection, "New York", 10, L2, bold, nohscroll, whiteBox;
	fontset = text, "New York", 12, L2, nohscroll;
	fontset = smalltext, "New York", 10, L2, nohscroll;
	fontset = input, "Courier", 12, L2, bold, nowordwrap;
	fontset = output, "Courier", 12, L2, nowordwrap;
	fontset = message, "Courier", 12, L2, R65535, nowordwrap;
	fontset = print, "Courier", 12, L2, nowordwrap;
	fontset = info, "Courier", 12, L2, nowordwrap;
	fontset = postscript, "Courier", 12, L2, nowordwrap;
	fontset = name, "Geneva", 10, L2, italic, B65535, nowordwrap, nohscroll;
	fontset = header, "Times", 10, L2;
	fontset = footer, "Times", 12, L2, center;
	fontset = help, "Geneva", 10, L2, nohscroll;
	fontset = clipboard, "New York", 12, L2;
	fontset = completions, "New York", 12, L2, nowordwrap;
	fontset = network, "Courier", 10, L2, nowordwrap;
	fontset = graphlabel, "Courier", 12, L2, nowordwrap;
	fontset = special1, "New York", 12, L2, nowordwrap;
	fontset = special2, "New York", 12, L2, center, nowordwrap;
	fontset = special3, "New York", 12, L2, right, nowordwrap;
	fontset = special4, "New York", 12, L2, nowordwrap;
	fontset = special5, "New York", 12, L2, nowordwrap;]
:[font = subsubtitle; inactive; ]
FortranDef by P.J. 14.06.90.
:[font = section; inactive; startGroup; ]
Usage
:[font = text; inactive; ]
 This file is PD software by
       Pekka Janhunen
       Finnish Meteorological Institute
       Geophysics Dept.
       pjanhune@finsun.csc.fi

Correct working at all conditions is not guaranteed!  
:[font = text; inactive; ]
This Mathematica code defines a function FortranDef which 
takes a Mathematica expression, writes an corresponding optimized
Fortran program, compiles it and defines a Mathematica function
which calls that program. Requires UNIX. Porting to other operating
systems (but not to Macintosh, of course) should be easy.
;[s]
7:0,0;5,1;16,0;67,1;78,0;168,1;179,0;315,-1;
3:4,17,12,Chicago,0,12,0,0,0;3,17,12,Chicago,2,12,0,0,0;0,17,12,Chicago,0,12,0,0,0;
:[font = input; initialization; ]
*)
BeginPackage["FortranDef`"]
(*
:[font = input; initialization; ]
*)
FortranDef::usage =
  "FortranDef[fn,expr,parlist,{opts}] takes expr, creates a corresponding Fortran program,
   compiles it and defines fn[parlist] so that it calls that program. Use VariableType -> 
   \"complex\" if your function contains complex numbers.";

FortranDef::notes =
  "System dependence: Modify the function CompileFortran if necessary."

MultiEval::usage =
 "MultiEval[f,{arglist}] evaluates f at many points. f must
  have been defined with FortranDef earlier.";

FortranIntegrate2DDef::usage =
  "FortranIntegrate2DDef[fn,f,{intvars},{othervars}] defines a
   two dimensional numerical integral routine in Fortran.\n
   Example: After the call\n\n
   \tFortranIntegrate2DDef[fn,f,{x,y},{u,v,..}]\n\n
   where f depends on x,y,u,v and w, the call\n\n
   \tFortranIntegrate2D[fn,{-1,1},{-1,1},{u0,v0,..}]\n\n
   gives the integral of f over the unit square and at 
   fixed u,v,..";

FortranIntegrate2D::usage =
  "Use FortranIntegrat2DDef prior to FortranIntegrate2D.\n
   Do ?FortranIntegrate2DDef.";

VariableType::usage =
  "VariableType is an option for FortranDef and FortranIntegrate2DDef.\n
   Use e.g. VariableType -> \"complex\".\n
   It determines the type of all variables used in internal computations.";

ParameterType::usage =
  "ParameterType is an option for FortranDef and FortranIntegrate2DDef.\n
   Use e.g. ParameterType -> "real*8".\n
   It determines the type of external parameters.";

IntvarType::usage =
  "IntvarType is an option for FortranDef and FortranIntegrate2DDef.\n
   Use e.g. IntvarType -> "real*4".\n
   It determines the type of integration variables.\n
   Making them complex makes no sense.";

Nx::usage =
  "Nx is an option for FortranIntegrate2DDef.\n
   Use e.g. Nx -> 100.\n
   It determines the number of grid points in x-direction.";

Ny::usage =
  "Ny is an option for FortranIntegrate2DDef.\n
   Use e.g. Ny -> 100.\n
   It determines the number of grid points in y-direction.";

DontCompile::usage =
  "DontCompile is an option for FortranDef and FortranIntegrate2DDef.\n
   Set DontCompile -> True to prevent unwanted recompilation\n
   and so to save time.";
(*
:[font = input; initialization; endGroup; ]
*)
Begin["`Private`"]
(*
:[font = section; inactive; startGroup; ]
Implementation
:[font = text; inactive; startGroup; whiteBox; ]
Define function FastSequence to give an "optimized" computation sequence
utilizing auxiliary variables...
:[font = input; initialization; endGroup; ]
*)
fast[a_?AtomQ]:= a;

FastOper[var_,expr_]:=
  If[ MemberQ[$rules,mySet[_,expr]],
      $rules[[ Position[$rules,mySet[_,expr]][[1,1]], 1 ]],
      AppendTo[$rules,mySet[var,expr]]; var ];
      
fast[u_ + v_]:= FastOper[Global`azxq[$Nrules+=1],fast[u]+fast[v]];
fast[u_ * v_]:= FastOper[Global`azxq[$Nrules+=1],fast[u]*fast[v]];
fast[u_ ^ v_]:= FastOper[Global`azxq[$Nrules+=1],fast[u]^fast[v]];
fast[f_[u_]]:= FastOper[Global`azxq[$Nrules+=1],f[fast[u]]]; (*general rule; use with caution !*)

FastSequence[f_]:=
  ($rules = {};
   $Nrules = 0;
   fast[f];
   {$rules,$Nrules});
(*
:[font = text; inactive; startGroup; whiteBox; ]
Modify FortranForm to work with list arguments. Enhance FortranForm
by giving rules for sqrt,power,asin etc.
:[font = input; initialization; endGroup; ]
*)
Unprotect[FortranForm];
FortranForm[x_List]:=
  OutputForm[SequenceForm @@ Flatten[{{FortranForm[#],","}& /@ Drop[x,-1], FortranForm[Last[x]]}]];
FortranForm[E^x_]:=
  OutputForm[SequenceForm["exp(",FortranForm[x],")"]];
FortranForm[x_^(1/2)]:=
  OutputForm[SequenceForm["sqrt(",FortranForm[x],")"]];
FortranForm[a_^b_]:=
  OutputForm[SequenceForm["(",FortranForm[a],")**(",FortranForm[b],")"]];
FortranForm[ArcSin[x_]]:=
  OutputForm[SequenceForm["asin(",FortranForm[x],")"]];
FortranForm[ArcCos[x_]]:=
  OutputForm[SequenceForm["acos(",FortranForm[x],")"]];
FortranForm[ArcTan[x_]]:=
  OutputForm[SequenceForm["atan(",FortranForm[x],")"]];
Protect[FortranForm];
(*
:[font = text; inactive; startGroup; whiteBox; ]
Define CompileFortran which calls the compiler. Define "take" which is the same
as Take but does not give error message if range is exceeded. Define also
FortranCall1, FortranCall and RatToReal that are needed later.
:[font = input; initialization; endGroup; ]
*)

CompileFortran[file_String]:=
  (Get[StringJoin["!f77 -o ",file,".prg ",file,".f"]]; Null);

take[list_List,n_Integer?Positive]:=
  If[n < Length[list], Take[list,n],list];

FortranCall[ProgramFileName_String,arglist_List,OutputType_String]:=
  Block[ {tmpfile,result,i},
    tmpfile = OpenTemporary[];
    Write[tmpfile,Length[arglist]];
    Write[tmpfile,FortranForm[#]]& /@ arglist;
    Close[tmpfile];
    result =
      ReadList[ StringJoin["!",ProgramFileName," <",tmpfile,
                           " | tr '(' ' ' | tr ')' ' ' | tr ',' ' '"],
                Number ];
    Get[StringJoin["!rm ",tmpfile]];
    If[StringMatchQ[StringJoin @@ (ToString /@ take[Characters[OutputType],7]),
                    "complex"],
       result = Table[result[[i]]+I result[[i+1]],
                      {i,1,Length[result],2}]];
    result]/;
  Count[Flatten[arglist],_?NumberQ] == Length[Flatten[arglist]];

FortranCall1[ProgramFileName_String,args_List,OutputType_String]:=
  (FortranCall[ProgramFileName,{args},OutputType] [[1]]) /;
  Count[args,_?NumberQ] == Length[args] ;

RatToReal[f_]:=
  (f//. Rational[n_,m_] :> N[n]/N[m]);
(*
:[font = text; inactive; startGroup; whiteBox; ]
Define FortranDef which writes the Fortran program and compiles it.
After the call FortranDef[fn,f,{x,y}] 
fn[x,y] gives the value of f at particular x,y.
Likewise, MultiEval[fn,{{x1,y1},{x2,y2},...}] yields a list where f is
evaluated at many points {xi,yi}.
:[font = input; initialization; endGroup; ]
*)
Options[FortranDef] =
  {VariableType -> "Real*8",
   ParameterType -> "Real*8",
   DontCompile -> False};

FortranDef[fn_,f_,vars_List,opts___]:=
  Block[
     {rules,Nrules,file,fileF,parlist,
      ProgramFileName,InputParams,OutputType,VariableType,ParameterType},
     VariableType = VariableType/.{opts}/.Options[FortranDef];
     ParameterType = ParameterType/.{opts}/.Options[FortranDef];
     parlist = Pattern[#,Blank[]]& /@ vars;
     file = ToString[fn];
     fileF = StringJoin[file,".f"];
     Clear[fn];
     OutputType = VariableType/.{opts}/.Options[FortranDef];
     If[!TrueQ[DontCompile/.{opts}/.Options[FortranDef]],
         {rules,Nrules} = FastSequence[f];
         OpenWrite[fileF];
         Write[fileF,OutputForm["\tprogram "], fn];
         Write[fileF,OutputForm["C Creator: Mathematica/FortranDef."]];
         Write[fileF,OutputForm["C Computes values of "],fn,
                     OutputForm[" in many points."] ];
         Write[fileF,OutputForm["\tinteger ij22f4"]];
         Write[fileF,OutputForm["\tinteger kk893g"]];
         Write[fileF,OutputForm[StringJoin["\t",VariableType," "]],
                     FortranForm[Global`azxq] ];
         Write[fileF,OutputForm[StringJoin["\t",ParameterType," "]],
                     FortranForm[vars] ];
         Write[fileF,OutputForm["\tdimension azxq("],Nrules+5,OutputForm[")"] ];
         Write[fileF,OutputForm["\tread(*,*) kk893g"]];
         Write[fileF,OutputForm["\tdo 10 ij22f4=1,kk893g"]];
         Write[fileF,OutputForm["\tread(*,*) "], FortranForm[vars] ];
         Write[fileF,OutputForm["\t"],
               FortranForm[#[[1]]], OutputForm["="],
               FortranForm[RatToReal[#[[2]]]]]& /@ rules;
         Write[fileF,OutputForm["\twrite(*,*) "],
                     FortranForm[Global`azxq[1]]];
         Write[fileF,OutputForm["10\tcontinue"]];
         Write[fileF,OutputForm["\tend"]];
         Close[fileF];
         CompileFortran[file] ];
     ProgramFileName = StringJoin["./",file,".prg"];
     InputParams = vars;
     arglist = .;
     MultiEval[fn,arglist_List] =
       FortranCall[ProgramFileName,arglist,OutputType];
     Release[fn @@ parlist] =
       FortranCall1[ProgramFileName, InputParams,OutputType];];
(*
;[s]
3:0,0;804,1;815,0;2232,-1;
3:2,14,10,Courier,1,12,0,0,0;1,12,10,Courier,3,12,0,0,0;0,12,10,Courier,0,12,0,0,0;
:[font = text; inactive; whiteBox; ]
Define FortranIntegrate2DDef.
After the call FortranIntegrate2DDef[fn,f,{x,y},{u,v,w}] 
FortranIntegrate2D[fn,{x1,x2},{y1,y2},{u0,v0,w0}] gives the definite
integral of f when x runs from x1 to x2, y runs from y1 to y2 and
the variables u,v,w have the particular numerical values u0,v0,w0.

FortranIntegrate2DDef writes a Fortran program that reads in 4+N+2
numbers (N=3 here) x1,x2; y1,y2; u0,v0,w0; Nx,Ny in this order, performs
the integral and writes the resulting single number to output. Nx and Ny
are the numbers of grid points in x and y directions, respectively. Semicolons
above indicate newlines.
:[font = input; initialization; ]
*)
Options[FortranIntegrate2DDef] =
  {VariableType -> "real*8",
   ParameterType -> "real*8",
   IntvarType -> "real*8",
   DontCompile -> False};

Options[FortranIntegrate2D] = {Nx -> 30, Ny -> 30};

FortranIntegrate2DDef\
  [fn_,integrand_,intvars_List,othervars_List,opts___]:=
  Block[
     {rules,Nrules,file,fileF,vars,opts2,
      ProgramFileName,InputParams,OutputType,
      variableType,parameterType,intvarType,
      x1,x2,y1,y2,others},
     variableType = VariableType/.{opts}/.Options[FortranIntegrate2DDef];
     parameterType = ParameterType/.{opts}/.Options[FortranIntegrate2DDef];
     intvarType = IntvarType/.{opts}/.Options[FortranIntegrate2DDef];
     parlist = If[TrueQ[Head[#]==Pattern],#,Pattern[#,Blank[]]]& /@ parlist;
     vars = #[[1]]& /@ parlist;
     file = StringJoin[ToString[fn],"Integrate2D"];
     fileF = StringJoin[file,".f"];
     OutputType = VariableType/.{opts}/.Options[FortranDef];
     If[!TrueQ[DontCompile/.{opts}/.Options[FortranDef]],
         {rules,Nrules} = FastSequence[integrand];
         OpenWrite[fileF];
         Write[fileF,OutputForm["\tprogram "], fn,OutputForm["INT2D"]];
         Write[fileF,OutputForm["C Creator: Mathematica/FortranIntegrate2DDef."]];
         Write[fileF,OutputForm["C Computes the 2D definite integral of "],fn];
         Write[fileF,OutputForm[StringJoin["\t",variableType," "]],
                     FortranForm[Global`azxq]];
         Write[fileF,OutputForm[StringJoin["\t",parameterType," "]],
                     FortranForm[othervars]];
         Write[fileF,OutputForm[StringJoin["\t",intvarType," "]],
                     FortranForm[intvars]];
         Write[fileF,
               OutputForm[StringJoin["\t",intvarType,
                                     " x1azxq,x2azxq,y1azxq,y2azxq,dxazxq,dyazxq"]]];
         Write[fileF,OutputForm[StringJoin["\t",variableType," sazxq"]]];
         Write[fileF,OutputForm["\tinteger Nxazxq,Nyazxq,Iazxq,Jazxq"]];
         Write[fileF,OutputForm["\tdimension "],
                     FortranForm[Global`azxq],
                     OutputForm["("],
                     Nrules+5,
                     OutputForm[")"] ];
         Write[fileF,OutputForm["C Read in integration limits, external parameters"]];
         Write[fileF,OutputForm["C and grid dimensions"]];
         Write[fileF,OutputForm["\tread(*,*) x1azxq,x2azxq"]];
         Write[fileF,OutputForm["\tread(*,*) y1azxq,y2azxq"]];
         Write[fileF,OutputForm["\tread(*,*) "], FortranForm[othervars]];
         Write[fileF,OutputForm["\tread(*,*) Nxazxq,Nyazxq"]];
         Write[fileF,OutputForm["\tsazxq=0."]];
         Write[fileF,OutputForm["\tdxazxq=(x2azxq-x1azxq)/Nxazxq"]];
         Write[fileF,OutputForm["\tdyazxq=(y2azxq-y1azxq)/Nyazxq"]];
         Write[fileF,OutputForm["C Start of integration loop"]];
         Write[fileF,OutputForm["\tdo 10 Iazxq=1,Nxazxq"]];
         Write[fileF,OutputForm["\tdo 20 Jazxq=1,Nyazxq"]];
         Write[fileF,OutputForm["\t"],
                     OutputForm[intvars[[1]]],
                     OutputForm["=x1azxq+dxazxq*(Iazxq-1)"]];
         Write[fileF,OutputForm["\t"],
                     OutputForm[intvars[[2]]],
                     OutputForm["=y1azxq+dyazxq*(Iazxq-1)"]];
         Write[fileF,OutputForm["\t"],
               FortranForm[#[[1]]], OutputForm["="],
               FortranForm[RatToReal[#[[2]]]]]& /@ rules;
         Write[fileF,OutputForm["\tsazxq=sazxq+"],
                     FortranForm[Global`azxq[1]]];
         Write[fileF,OutputForm["20\tcontinue"]];
         Write[fileF,OutputForm["10\tcontinue"]];
         Write[fileF,OutputForm["\twrite(*,*) sazxq*dxazxq*dyazxq"]];
         Write[fileF,OutputForm["\tend"]];
         Close[fileF];
         CompileFortran[file] ];
     ProgramFileName = StringJoin["./",file,".prg"];
     FortranIntegrate2D
       [fn,{x1_Real,x2_Real},{y1_Real,y2_Real},
        others_List?(Length[#]==Length[othervars]&),opts2___] =
         IntegrateCall[ProgramFileName,x1,x2,y1,y2,others,opts2];];

IntegrateCall\
  [ProgramFileName_String,
   x1_Real,x2_Real,y1_Real,y2_Real,others_List,opts___]:=
  Block[ {tmpfile,result,i,nx,ny},
    nx = Nx/.{opts}/.Options[FortranIntegrate2D];
    ny = Ny/.{opts}/.Options[FortranIntegrate2D];
    tmpfile = OpenTemporary[];
    Write[tmpfile,FortranForm[{x1,x2}]];
    Write[tmpfile,FortranForm[{y1,y2}]];
    Write[tmpfile,FortranForm[others]];
    Write[tmpfile,FortranForm[{nx,ny}]];
    Close[tmpfile];
    result =
      ReadList[ StringJoin["!",ProgramFileName," <",tmpfile,
                           " | tr '(' ' ' | tr ')' ' ' | tr ',' ' '"],
                Number ];
    Get[StringJoin["!rm ",tmpfile]];
    If[StringMatchQ[StringJoin @@ (ToString /@ take[Characters[OutputType],7]),
                    "complex"],
       result = Table[result[[i]]+I result[[i+1]],
                      {i,1,Length[result],2}]];
    result[[1]]]/;
  Count[Flatten[others],_?NumberQ] == Length[Flatten[others]];

(*
;[s]
4:0,0;1178,1;1189,0;4106,2;4958,-1;
3:2,14,10,Courier,1,12,0,0,0;1,12,10,Courier,3,12,0,0,0;1,14,10,Courier,1,12,0,0,0;
:[font = text; inactive; startGroup; whiteBox; ]
End package
:[font = input; initialization; ]
*)
Print[{FortranDef,MultiEval,FortranIntegrate2DDef,FortranIntegrate2D}]
(*
:[font = input; initialization; ]
*)
End[]
(*
:[font = input; initialization; endGroup; endGroup; ]
*)
EndPackage[]
(*
:[font = section; inactive; startGroup; Cclosed; ]
Example
:[font = input; startGroup; ]
expr = D[ Exp[Sqrt[x^2+y^4]] Sin[x^2 + y^4], {x,0},{y,0}];
Clear[FortranExpr];
Expr[x_,y_] = expr;
FortranDef[FortranExpr, expr, {x_,y_}, VariableType -> "real*8" ]
:[font = print; inactive; endGroup; endGroup; ]

:[font = section; inactive; startGroup; Cclosed; ]
Testing
:[font = input; ]
inp = Table[{x},{x,-1.,1.,.1}];
:[font = input; startGroup; ]
Timing[res = (Expr /@ inp);]
:[font = output; inactive; output; endGroup; ]
{6.700000000000001*Second, Null}
;[o]
{6.7 Second, Null}
:[font = input; startGroup; ]
Timing[fres = MultiEval[FortranExpr,inp];]
:[font = output; inactive; output; endGroup; ]
{0.43*Second, Null}
;[o]
{0.43 Second, Null}
:[font = input; ]
Table[(fres[[i]]-res[[i]])/res[[i]],{i,1,Length[res]}]
:[font = input; ]
res
:[font = input; endGroup; ]
fres
^*)