How to change an expression to Fortran90 form?

Question

I have a large expression produced by Mathematica. To use this expression in my Fortran program I need to change it to a Fortran90 form. I am not satisfied with the FortranForm command.

The output of the Mathematica expression is like this,

       -48*Sqrt[6]*E^((-2*S[1, x]^2)/bx - (S[1, y]^2)/bx - (S[1, z]^2)/
        bz - (S[2, x]^2)/bx - (S[2, y]^2)/bx - (S[2, z]^2)/
        bz - (Sp[1, x]^2)/bx - (Sp[1, y]^2)/bx - (Sp[1, z]^2)/bz)
         + 3*Sqrt[5*v]

I hope I can get this expression like this,

    -48.q0*sqrt(6.q0)*Exp((-(2.0*s1x**2/bx)-s1y**2/bx-s1z**2/bz  &
    -s2x**2/bx-s2y**2/bx - s2z**2/bz - sp1x**2/bx - sp1y**2/bx   &
    - sp1z**2/bz)+ 3.q0*Sqrt[5.q0*v]

Notice that:

1. The width should be fixed and there is a "&" at the end of the line. It is better the single line space can be fully occupied.

2. Some integers should be changed to reals. For example, 6 -> 6.q0 or 6.d0;

3. It is better for E^ to be changed to Exp() form.

How can I change S[a,b] to Sab easily? This is not a big problem. I could replace it in Mathematica.

Anyway, I want to change my Mathematica expression into a Fortran90 expression automaticly. Is it possible?

Answer 1

For simplicity and elegance I use Verbeia's nice Riffle-ing for continuation. While it is possible to format real numbers as requested in FortranForm I always thought that this is way too complicated in Mathematica (the code below is slightly modified from my old FeynCalc-Write2 function described here

myFortrantoggle=False; (* needed to avoid recursion *)
Clear[q];
Unprotect[Real];
Real/:Format[r_Real/;r>=0,FortranForm]:=({mantissa,exponent}=MantissaExponent[r];
If[r===0.,exponent=1];
SequenceForm[10. mantissa,q,exponent-1])/;(myFortrantoggle=!myFortrantoggle);
Protect[Real];
Unprotect[Times];
ClearAttributes[Times,Orderless];(* not really needed, but keeps order *)
step1=FortranForm[(-48*Sqrt[6]*E^((-2*S[1,x]^2)/bx-(S[1,y]^2)/bx-(S[1,z]^2)/bz-(S[2,x]^2)/bx-(S[2,y]^2)/bx-(S[2,z]^2)/bz-(Sp[1,x]^2)/bx-(Sp[1,y]^2)/bx-(Sp[1,z]^2)/bz)+3*Sqrt[5*v])/.
{Sqrt[a_]:>sqrt[N@a],Sqrt[i_Integer]Sqrt[f_Symbol]:>sqrt[N@i f],E^b_ :> exp[b]}/.
(j_Integer*x_):>(N[j]x)/.
{S[a_,b_]:>ToExpression[StringJoin["s",ToString[a],ToString[b]]],Sp[a_,b_]:>ToExpression[StringJoin["sp",ToString[a],ToString[b]]]}];
result=StringJoin@@Riffle[With[{splits=StringSplit[ToString@step1," "]},Fold[If[StringLength[Last@#1]+StringLength[#2]>60,Join[#1,{#2}],Join[Most[#1],{StringJoin[Last[#1],#2]}]]&,{First@splits},Rest[splits]]],"  &\n"];
SetAttributes[Times,Orderless];Protect[Times];
result

Answer 2

You can get a long way to your goal using FortranForm in Mathematica and some replacement rules. ToExpression and ToString are very useful for this kind of application. (I'm giving this one a name for later use below.)

step1 = FortranForm[(-48*Sqrt[6]*
  E^((-2*S[1, x]^2)/bx - (S[1, y]^2)/bx - (S[1, z]^2)/
     bz - (S[2, x]^2)/bx - (S[2, y]^2)/bx - (S[2, z]^2)/
     bz - (Sp[1, x]^2)/bx - (Sp[1, y]^2)/bx - (Sp[1, z]^2)/bz) + 
  3*Sqrt[5*v])] /. {S[a_, b_] :> ToExpression[StringJoin["s", ToString[a], ToString[b]]], 
  Sp[a_, b_] :>  ToExpression[StringJoin["sp", ToString[a], ToString[b]]]}

-48*Sqrt(6)*E**
 -    ((-2*s1x**2)/bx - s1y**2/bx - s1z**2/bz - s2x**2/bx - s2y**2/bx - 
 -      s2z**2/bz - sp1x**2/bx - sp1y**2/bx - sp1z**2/bz) + 
 -  3*Sqrt(5)*Sqrt(v)

However I don't think you can get your preferred format for integers using Mathematica alone. If you try to append a period (dot) to an integer expression, Mathematica will try to calculate the actual real-valued expression as shown, treating q0 as a separate symbolic quantity:

FortranForm[(-48*Sqrt[6]*
     E^((-2*S[1, x]^2)/bx - (S[1, y]^2)/bx - (S[1, z]^2)/
         bz - (S[2, x]^2)/bx - (S[2, y]^2)/bx - (S[2, z]^2)/
         bz - (Sp[1, x]^2)/bx - (Sp[1, y]^2)/bx - (Sp[1, z]^2)/bz) + 
  3*Sqrt[5*v])] /. {S[a_, b_] :> ToExpression[StringJoin["s", ToString[a], ToString[b]]], 
  Sp[a_, b_] :> ToExpression[StringJoin["sp", ToString[a], ToString[b]]],
  z_Integer :> ToExpression[ToString[z]~StringJoin~".q0"]} 

   -117.57550765359254*E**
 -    ((-2.*q0*s1x**(2.*q0))/bx**(1.*q0) - 
 -      (1.*q0*s1y**(2.*q0))/bx**(1.*q0) - 
 -      (1.*q0*s1z**(2.*q0))/bz**(1.*q0) - 
 -      (1.*q0*s2x**(2.*q0))/bx**(1.*q0) - 
 -      (1.*q0*s2y**(2.*q0))/bx**(1.*q0) - 
 -      (1.*q0*s2z**(2.*q0))/bz**(1.*q0) - 
 -      (1.*q0*sp1x**(2.*q0))/bx**(1.*q0) - 
 -      (1.*q0*sp1y**(2.*q0))/bx**(1.*q0) - 
 -      (1.*q0*sp1z**(2.*q0))/bz**(1.*q0))*q0**1.5 + 
 -  6.708203932499369*q0**1.5*Sqrt(v)

Similarly, I'm not sure there is a programmatic way to convert the E** form to exp(...) form without building a more complex parser.

Here is a partial solution to your need for ampersands at the end of the line. The basic idea is to convert to a string, split at the appropriate length, and then Riffle in the ampersand and a line break. However, you can't just break at the 80-character limit because that breaks expressions like sp1z across lines, which is not what you want. Probably what you should try is to split the string progressively at whitespace, and then build up piece-by-piece into strings of up to (but no more than) the desired number characters using Fold, like this (I've used 60 characters for illustrative purposes):

StringJoin @@ 
 Riffle[With[{splits = StringSplit[ToString@step1, " "]}, 
   Fold[If[StringLength[Last@#1] + StringLength[#2] > 60, 
      Join[#1, {#2}], 
      Join[Most[#1], {StringJoin[Last[#1], #2]}]] &, {First@splits}, 
    Rest[splits]] ], "  &\n"]

-48*Sqrt(6)*E**((-2*s1x**2)/bx-s1y**2/bx-s1z**2/bz-s2x**2/bx  &
  -s2y**2/bx-s2z**2/bz-sp1x**2/bx-sp1y**2/bx-sp1z**2/bz)+  &
  3*Sqrt(5)*Sqrt(v)

Answer 3

This changes E**argument to exp(argument):

Unprotect[Power];
Power /: Format[Power[E, x_], FortranForm] := exp[x]
Protect[Power];
{FortranForm[E^(3*z)], FortranForm[Exp[2*z^3]]}

(* Out: {exp(3*z), exp(2*z**3)} *)