(* This file demonstrates the use of the FortranDef -package.
   The package assumes that your Mathematica is running under
   UNIX. Modifications to other operating systems (but not to
   Macintosh, of course) should be easy.
   
   FortranDef.m is PD software by
             Pekka Janhunen
             Finnish Meteorological Institute
             Geophysics Department
             pjanhune@finsun.csc.fi  *)

(* First we load in the package..*)
<<Misc/FortranDef.m

?FortranDef

(* Let's define a complicated function f:*)
f[k0_,x_,y_] = D[ Exp[-I k0 Sqrt[x^2+y^2]],{x,4},{y,4}];

(* Then the same function with Fortran:*)
Clear[Ff];
FortranDef[Ff,f[k0,x,y],{k0_,x_,y_},VariableType->"complex*16"];

f[0.23,54.,34.]  (* this is the Mathematica function *)

Ff[0.23,54.,34.]  (* and this calls the Fortran version *)

(* Some test data *)
inps = Table[{k0,23.54,34.},{k0,0.1,0.9,0.01}];

(* This calls the Fortran program for each entry:*)
Timing[ Ffres = ((Ff @@ #)& /@ inps); ]

(* And this calls the Fortran code only once for the whole list:*)
Timing[ Fres = MultiEval[Ff,inps]; ]

(* For comparison, we do the same with Mathematica version: *)
Timing[ res = ((f @@ #)& /@ inps); ]

(* Verify MultiEval result and Mathematica result:*)
Max[Table[ Abs[(res[[i]] - Fres[[i]])/res[[i]]], {i,1,Length[inps]}]]

(* Verify Ff result and Mathematica result:*)
Max[Table[ Abs[(res[[i]] - Ffres[[i]])/res[[i]]], {i,1,Length[inps]}]]

(* We can do any Mathematica operations with Ff: *)
Plot[Re[Ff[k0,10.3,40.3]],{k0,0.1,0.9}]

(* Plotting etc. becomes even faster if we use MultiEval
   to compute the values. *)
inps = Table[{k0,10.3,40.3},{k0,0.1,0.9,0.01}];
vals = MultiEval[Ff,inps];
ListPlot[Re[vals], PlotJoined -> True]

(* ByteCount tells us how large the Mathematica expression for
   f actually is:*)
 ByteCount[ f[k0,x,y] ]

(* The generated Fortran program listing is also large: *)
Run[ "cat Ff.f | wc" ]