How to change an expression to Fortran90 form?
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
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)
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)} *)