(*^

::[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. 25.05.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!  

 Note: global variables:
    - $rules
    - $Nrules
    - $OutputType
:[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; endGroup; ]
*)
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.";
(*
:[font = section; inactive; startGroup; ]
Definitions
:[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[a[$Nrules+=1],fast[u]+fast[v]];
fast[u_ * v_]:= FastOper[a[$Nrules+=1],fast[u]*fast[v]];
fast[u_ ^ v_]:= FastOper[a[$Nrules+=1],fast[u]^fast[v]];
fast[f_[u_]]:= FastOper[a[$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. Define CompileFortran and FortranDef...
:[font = input; initialization; endGroup; 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];

CompileFortran[file_String]:=
  Run["f77", "-o", StringJoin[file,".prg"], StringJoin[file,".f"]];

FortranCall[ProgramFileName_String,arglist_List]:=
  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 ];
    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]:=
  (FortranCall[ProgramFileName,{args}] [[1]]) /;
  Count[args,_?NumberQ] == Length[args] ;

RatToReal[f_]:=
  (f//. Rational[n_,m_] :> N[n]/N[m]);

Options[FortranDef] = {VariableType -> "Real*8", ParameterType -> "Real*8"};

FortranDef[fn_,f_,parlist_List,opts___]:=
  Block[ {rules,Nrules,file,fileF,vars,
          ProgramFileName,InputParams},
         {rules,Nrules} = FastSequence[f];
         vars = #[[1]]& /@ parlist;
         file = ToString[fn];
         fileF = StringJoin[file,".f"];
         OpenWrite[fileF];
         $OutputType = VariableType/.opts/.Options[FortranDef];
         Write[fileF,OutputForm["\tprogram "], fn];
         Write[fileF,OutputForm["\tinteger ij22f4"]];
         Write[fileF,OutputForm["\tinteger kk893g"]];
         Write[fileF,OutputForm[StringJoin["\t",VariableType/.opts/.Options[FortranDef]," "]],a];
         Write[fileF,OutputForm[StringJoin["\t",ParameterType/.opts/.Options[FortranDef]," "]],FortranForm[vars]];
         Write[fileF,OutputForm["\tdimension a("],Nrules+5,OutputForm[")"] ];
         Write[fileF,OutputForm["\tread(*,*) kk893g"]];
         Write[fileF,OutputForm["\tdo 111 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[a[1]]];
         Write[fileF,OutputForm["111\tcontinue"]];
         Write[fileF,OutputForm["\tend"]];
         Close[fileF];
         CompileFortran[file];
         ProgramFileName = StringJoin["./",file,".prg"];
         InputParams = vars;
         arglist = .;
         MultiEval[fn,arglist_List] = FortranCall[ProgramFileName,arglist];
         Release[fn @@ parlist] = FortranCall1[ProgramFileName, InputParams];];
(*
:[font = section; inactive; startGroup; ]
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; ]
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
^*)