5
$\begingroup$

I need to transform a function from Mathematica to Fortran77. An internal program FortranForm gives some result, but it works incorrectly in Fortran77. The string must start a 7th column and finish at 72th, ArcTan and Log must be transformed to ATAN and DLOG respectively, numbers also must be transformed (e.g. 2. -> 2d0).

I have a function:

 C3 (c - y) + 
 C3 (c + y) + ((1 + nu) (-1 + 2 nu) P (c - y) ArcTan[(-c + y)/
   x])/(π Y) + ((1 + nu) (-1 + 2 nu) P (c + y) ArcTan[(c + y)/
   x])/(π Y) + ((-1 + nu) (1 + nu) P x Log[
   x^2 + (c - y)^2])/(π Y) + ((1 + nu) P (x - nu x) Log[
   x^2 + (c + y)^2])/(π Y)

I'd like to have the following result:

< V = C3*(c - y) + C3*(c + y) + 
 &  ((1D0 + nu)*(-1D0 + 2D0*nu)*P*(c - y)*ATAN((-c + y)/x))/(Pi*E) + 
 &  ((1D0 + nu)*(-1D0 + 2D0*nu)*P*(c + y)*ATAN((c + y)/x))/(Pi*E) + 
 &  ((-1D0 + nu)*(1D0 + nu)*P*x*DLOG(x**2 + (c - y)**2))/(Pi*E) + 
 &  ((1D0 + nu)*P*(x - nu*x)*DLOG(x**2 + (c + y)**2))/(Pi*E)
>       

Here I have first string starting from 7-th columns (V=...), symbol & at 6-th position show continuation, 1->1D0, ArcTan->ATAN, Y->E, etc.

$\endgroup$
  • 1
    $\begingroup$ "a function" - where is it? $\endgroup$ – J. M. will be back soon Mar 13 '16 at 14:06
  • 5
    $\begingroup$ edit the question with the details. Why f77? Every fortran compiler in the past 20 years supports f90 at least. Use freeform and forget about the column 7 issue. $\endgroup$ – george2079 Mar 13 '16 at 15:45
  • 1
    $\begingroup$ FortranForm cannot translate functions. It can only translate simple expressions, and even then its result is more of a starting point ... I would use ReplaceAll on ArcTan -> ATAN and would edit the result with a text editor (including fixing numbers). $\endgroup$ – Szabolcs Mar 13 '16 at 16:55
  • 1
    $\begingroup$ Please add a useful example of input/output to your question to work with. The more specific your question, the higher the chance to get good answers. $\endgroup$ – Yves Klett Mar 13 '16 at 17:06
  • 3
    $\begingroup$ FORTRAN77 Code Generation Package library.wolfram.com/infocenter/Demos/60 $\endgroup$ – Orders Mar 14 '16 at 10:50
9
$\begingroup$

Warning: as @george points out in the comments, something is wrong with the output of this package. It produced incorrect d0d0. I don't have time to investigate at the moment, so making it community wiki. Feel free to post your own (corrected) answer based on the information here.

Using the Format.m package that @Orders posted,

expr = C3 (c - y) + 
   C3 (c + y) + ((1 + nu) (-1 + 2 nu) P (c - y) ArcTan[(-c + y)/
        x])/(π Y) + ((1 + nu) (-1 + 2 nu) P (c + 
        y) ArcTan[(c + y)/x])/(π Y) + ((-1 + nu) (1 + nu) P x Log[
       x^2 + (c - y)^2])/(π Y) + ((1 + nu) P (x - nu x) Log[
       x^2 + (c + y)^2])/(π Y);

FortranAssign["expr", expr, AssignReplace -> {"log" -> "dlog"}, AssignOptimize -> False]

gives the following

        expr = C3*(c - y) + C3*(c + y) + (3.183098861837907d-1*(1.d0 + n
     &  u)*(-1.d0d0 + 2.d0*nu)*P*(c - y)*atan((-c + y)/x))/Y + (3.183098
     &  861837907d-1*(1.d0 + nu)*(-1.d0d0 + 2.d0*nu)*P*(c + y)*atan((c +
     &   y)/x))/Y + (3.183098861837907d-1*(-1.d0d0 + nu)*(1.d0 + nu)*P*x
     &  *dlog((c - y)**2 + x*x))/Y + (3.183098861837907d-1*(1.d0 + nu)*P
     &  *(x - nu*x)*dlog((c + y)**2 + x*x))/Y

The 3.183d-1 you see here is simply a decimal approximation of 1/Pi to double precision.

We can also use automatic expression optimization, which changes the formula to an equivalent form that minimizes the number of arithmetic operations. For this it needs to use temporary variables, which I named tempXX here.

FortranAssign["expr", expr, AssignReplace -> {"log" -> "dlog"}, 
 "OptimizationSymbol" -> temp]

        temp1 = -y
        temp2 = c + temp1
        temp6 = 1.d0 + nu
        temp7 = 2.d0*nu
        temp8 = -1.d0d0 + temp7
        temp4 = c + y
        temp9 = 1/Y
        temp10 = 1/x
        temp20 = x*x
        expr = C3*temp2 + C3*temp4 + 3.183098861837907d-1*P*temp4*temp6*
     &  temp8*temp9*atan(temp10*temp4) + 3.183098861837907d-1*P*temp2*te
     &  mp6*temp8*temp9*atan(temp10*(-c + y)) + 3.183098861837907d-1*(-1
     &  .d0d0 + nu)*P*temp6*temp9*x*dlog(temp20 + temp2*temp2) + 3.18309
     &  8861837907d-1*P*temp6*temp9*(x - nu*x)*dlog(temp20 + temp4*temp4
     &  )

Getting the package to work in Mathematca 10:

If you try to load it as is, it will complain about not finding Utilities`FilterOptions`, a package that no longer ships with Mathematica 10. To fix this,

  • Change BeginPackage["Format`", "Utilities`FilterOptions`"] to BeginPackage["Format`"].

  • Add FilterOptions[fun_, opts___] := Sequence@@FilterRules[{opts}, Options[fun]] right after Begin["Private`"].

$\endgroup$
  • $\begingroup$ looks a little buggy, d0d0 ? $\endgroup$ – george2079 Mar 14 '16 at 14:57
  • $\begingroup$ I use Mathematica 10.3. I've got the package Format.m, changed it using instructions and installed. However when I try to use a command <&lt;Format.m from this example packagedata.net/index.php/links/examples/id/207 , I have a following mistake "Syntax::sntxb: Expression cannot begin with "<&lt". " What can be wrong? $\endgroup$ – IgorPy Mar 15 '16 at 8:44
  • $\begingroup$ @IgorPy That should be <<Format.m. However, as you can see, the package misformats negative numbers and gives -1.d0d0 with a double d0d0. I looked at how to fix this and I came away with the conclusion that it's not worth it. It depends on internal Mathematica implementation details that are messy and changed since this package was written. I think you're better off manually editing the output and fixing up problems. Unless you really need to translate hundreds and hundreds of expressions in which case you can write your own package. Give that at least 1-2 days of time. $\endgroup$ – Szabolcs Mar 15 '16 at 8:53
  • $\begingroup$ @Szabolcs, thanks a lot! I will use your package and change d0d0 to d0 manually. Now it works. $\endgroup$ – IgorPy Mar 15 '16 at 9:10
  • $\begingroup$ @IgorPy It's not my package ... if it were I would fix it ... $\endgroup$ – Szabolcs Mar 15 '16 at 9:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.