(*^

::[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", 10, L2, nohscroll;
	fontset = smalltext, "New York", 9, L2, nohscroll;
	fontset = input, "Courier", 10, L2, bold, nowordwrap;
	fontset = output, "Courier", 10, L2, nowordwrap;
	fontset = message, "Courier", 10, L2, R65535, nowordwrap;
	fontset = print, "Courier", 10, L2, nowordwrap;
	fontset = info, "Courier", 10, 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 = input; ]
(* 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
:[font = input; startGroup; ]
?FortranDef
:[font = info; inactive; endGroup; ]
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.
:[font = input; startGroup; ]
(* 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"];
:[font = print; inactive; endGroup; ]

:[font = input; startGroup; ]
f[0.23,54.,34.]  (* this is the Mathematica function *)
:[font = output; inactive; output; endGroup; ]
-(1.660300910198964*10^-7) - 3.087544567735867*10^-7*I
;[o]
          -7             -7
-1.6603 10   - 3.08754 10   I
:[font = input; startGroup; ]
Ff[0.23,54.,34.]  (* and this calls the Fortran version *)
:[font = output; inactive; output; endGroup; ]
-(1.6603009101989165*10^-7) - 3.0875445677358941*10^-7*I
;[o]
                      -7                        -7
-1.6603009101989165 10   - 3.0875445677358941 10   I
:[font = input; ]
(* Some test data *)
inps = Table[{k0,23.54,34.},{k0,0.1,0.9,0.01}];
:[font = input; startGroup; ]
(* This calls the Fortran program for each entry:*)
Timing[ Ffres = ((Ff @@ #)& /@ inps); ]
:[font = output; inactive; output; endGroup; ]
{7.280000000000001*Second, Null}
;[o]
{7.28 Second, Null}
:[font = input; startGroup; ]
(* And this calls the Fortran code only once for the whole list:*)
Timing[ Fres = MultiEval[Ff,inps]; ]
:[font = output; inactive; output; endGroup; ]
{3.39*Second, Null}
;[o]
{3.39 Second, Null}
:[font = input; startGroup; ]
(* For comparison, we do the same with Mathematica version: *)
Timing[ res = ((f @@ #)& /@ inps); ]
:[font = output; inactive; output; endGroup; ]
{52.29*Second, Null}
;[o]
{52.29 Second, Null}
:[font = input; startGroup; ]
(* Verify MultiEval result and Mathematica result:*)
Max[Table[ Abs[(res[[i]] - Fres[[i]])/res[[i]]], {i,1,Length[inps]}]]
:[font = output; inactive; output; endGroup; ]
1.14355601322054*10^-12
;[o]
          -12
1.14356 10
:[font = input; startGroup; ]
(* Verify Ff result and Mathematica result:*)
Max[Table[ Abs[(res[[i]] - Ffres[[i]])/res[[i]]], {i,1,Length[inps]}]]
:[font = output; inactive; output; endGroup; ]
1.14355601322054*10^-12
;[o]
          -12
1.14356 10
:[font = input; startGroup; ]
(* We can do any Mathematica operations with Ff: *)
Plot[Re[Ff[k0,10.3,40.3]],{k0,0.1,0.9}]
:[font = output; inactive; output; endGroup; ]
The Unformatted text for this cell was not generated.
Use options in the Actions Settings dialog box to control
when Unformatted text is generated.
;[o]
-Graphics-
:[font = input; startGroup; ]
(* 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]
:[font = output; inactive; output; endGroup; ]
The Unformatted text for this cell was not generated.
Use options in the Actions Settings dialog box to control
when Unformatted text is generated.
;[o]
-Graphics-
:[font = input; startGroup; ]
(* ByteCount tells us how large the Mathematica expression for
   f actually is:*)
 ByteCount[ f[k0,x,y] ]
:[font = output; inactive; output; endGroup; ]
40316
;[o]
40316
:[font = input; startGroup; ]
(* The generated Fortran program listing is also large: *)
Run[ "cat Ff.f | wc" ]
:[font = print; inactive; ]
       312       432      651
:[font = output; inactive; output; endGroup; ]
0
;[o]
0
^*)