Might be very naive; I am trying for the first time to develop a Mathematica package as shown below. The intended usage is that the user will load the package, specify the input functions 'Model' and 'ShapeFunc', then use the 'PlotWEC' function from the package to generate a plot for 'WEC' eveluated by the package. My package loads, but I get an erroneous plot intersecting the x axis, perhaps meaning that WEC is evaluated as 0 for some reason. Here is the package code:


PlotWEC::usage = 
        " PlotWEC[Model_,ShapeFunc_] Plots the behavior of the WEC with r for the given modified gravity model and shape function";
PlotWEC[Model_,ShapeFunc_]:=Module[{FDModel= D[Model,R],R=(2 D[ShapeFunc, r])/r^2,WEC=(FDModel*D[ShapeFunc,r])/r^2}, Plot [{WEC},{r,0.1,10},PlotTheme->"Scientific", PlotStyle->{Blue}, LabelStyle->{FontSize->12,FontFamily->"Arial", Black,Bold},FrameLabel->{"r","\[Rho]"},  ImageSize-> Medium]];



Also, here is my usage to test the package:

Import[NotebookDirectory[] <> "CWAP.wl"]

Model = 0.1*R^2;

ShapeFunc =  (0.9 Log[r + 1])/Log[0.9 + 1];

PlotWEC[Model, ShapeFunc]

The 'R' is contained in the model defined. FDModel is the first derivative of Model with respect to R. And when evaluating WEC, R is taken as a derivative of 'ShapeFunc' wrt 'r'. I think my problem is in the definition of R. Please suggest a working code for the package. It would be really helpful if someone can explain where I am wrong. Thanks in advance!

EDIT: Here is the guide notebook I made to develop the package.


(*Define the f(R) gravity model used. This will be specified by the user when using the package*)

fRModel = 0.1 R^2;

(*Define the desired shape function. This will be specified by the user when using the package*)

ShapeFunc = (q Log[r + 1])/Log[q + 1];

(*Define the Ricci Scalar*)

R = (2 D[ShapeFunc, r])/r^2;

(*Specify the location of the throart. This will be specified by the user when using the package*)

q = 0.9;

(*The expression for the Weak Energy Condition is *)

WEC = (D[fRModel,r] * D[ShapeFunc, r])/r^2;

(*Plot the WEC*)

GraphWEC = 
  Plot[WEC, {r, 0.1, 20}, PlotTheme -> "Scientific", 
   PlotStyle -> {Blue}, 
   LabelStyle -> {FontSize -> 12, FontFamily -> "Arial", Black, Bold},
    FrameLabel -> {"r", "\[Rho]"}, ImageSize -> Medium];

The user should be able to input the Model and the ShapeFunction and get the plot for the WEC. In my previous code for the package implementation, I had specified 'q' manually as 0.9. Please note that in evaluating WEC, once D[fRModel,r] has been evaluated, and remaining 'R's should be evaluated as R = (2 D[ShapeFunc, r])/r^2.


1 Answer 1


Passing "raw" symbols around is fraught. One thing that's happening here is that the symbol R used in PlotWEC is not the same symbol as the one you use to define Model. A new, unqualified (meaning it uses no explicit context) symbol will be assigned to the context that is currently in effect. When you defined PlotWEC, the current context was CWAP`Private` . When you defined Model, the current context was Global` . One thing you could try would be to use functions rather than "raw" expressions. You can use Derivative to take derivatives of functions.

Making this concrete...

  1. Change your definition of PlotWEC to take functions instead of expressions (this would still be between the Begin and End).

    PlotWEC[Model_, ShapeFunc_] :=
        (Derivative[1][Model][r]*Derivative[1][ShapeFunc][r])/r^2, {r, 0.1, 10},
        <...put your options here...>]
  2. Change your Model and ShapeFunc to be functions.

    Model[r_] := 0.1*r^2;
    ShapeFunc[r_] := (0.9 Log[r + 1])/Log[0.9 + 1];
  3. Call PlotWEC with these functions.

    PlotWEC[Model, ShapeFunc]

I didn't actually try to analyze and understand what you're doing, so I don't know if the resulting plot is what you were expecting. There may be other things to fix.


(FWIW, the above refactoring produces the wrong answer because I initially misunderstood the purpose of capital R.)

The following demonstrates how to eliminate the dependence on global variables. You ought to be able to apply this pattern in your package definition.

fRModel[R_] := 0.1 R^2;
F1 = Derivative[1][fRModel];
F2 = Derivative[2][fRModel];(* this doesn't seem to be used anywhere, but I've kept it for future reference *)
ShapeFunc[r_, q_] := (q Log[r + 1])/Log[q + 1];
R[r_, q_] := (2 Derivative[1, 0][ShapeFunc][r, q])/r^2;
WEC[r_, q_] := (F1[R[r, q]]*Derivative[1, 0][ShapeFunc][r, q])/r^2;
  WEC[r, 0.9], {r, 0.1, 20}, 
  PlotTheme -> "Scientific", 
  PlotStyle -> {Blue}, 
  LabelStyle -> {FontSize -> 12, FontFamily -> "Arial", Black, Bold}, 
  FrameLabel -> {"r", "\[Rho]"}, 
  ImageSize -> Medium]
  • $\begingroup$ Thank you very much for the answer. One small thing, 'r' and 'R' are different things; 'r' is the radial coordinate, 'R' is the Ricci scalar. Model[R_]:=R^2 would be the correct definition. I am trying to plot the behavior of the 'weak energy condition'. The correct plot should give some positive values going to zero asympotically (as verified before). I tried implementing your changes, but using PlotWEC[Model, ShapeFunc] now yields an empty plot. One more thing, when we evaluate Derivative[1][Model][R], we get 0.2 R; and therein, R should take the value D[ShapeFunc,r]/r^2 $\endgroup$
    – madmiKe
    Mar 8, 2023 at 17:52
  • $\begingroup$ There are probably some stray symbol definitions left around, so try killing the kernel and evaluating everything again. $\endgroup$
    – lericr
    Mar 8, 2023 at 17:54
  • $\begingroup$ As for r versus R, yeah, I hadn't grokked that before. I don't know anything about the context nor the vocabulary that you're using, so if you still need help you'll need to explicitly provide the exact formulas and definitions you want--I won't be able to just figure it out. $\endgroup$
    – lericr
    Mar 8, 2023 at 17:57
  • $\begingroup$ I am providing you the guide notebook I made to develop the package in the first place. Please see my edits to the question. I am sure you will understand the issue. $\endgroup$
    – madmiKe
    Mar 8, 2023 at 17:58
  • $\begingroup$ Please check the code at the end of the question. Thanks! $\endgroup$
    – madmiKe
    Mar 8, 2023 at 18:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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