# Difficulty plotting sum of functions

I'm interested in plotting this function (provided in equation 11 of this paper):

\begin{aligned} \begin{aligned} S_{S}(\omega)=& \sum_{m_{5}, m_{I}, m_{S}^{\prime}, m_{I}^{\prime}=-\infty}^{\infty} \operatorname{sinc}^{2}\left[\frac{\tau_{0}}{2}\left(\omega-\omega_{S}\right)\right] \\ & \times \frac{\left(\frac{\gamma_{S}}{2}\right)^{2}+\left(\omega-\left(\omega_{S}+m_{S}^{\prime} \Delta \omega_{S}\right)\right)\left(\omega-\left(\omega_{S}+m_{S} \Delta \omega_{S}\right)\right)}{\left(\left(\frac{\gamma_{S}}{2}\right)^{2}+\left(\omega-\left(\omega_{S}+m_{S}^{\prime} \Delta \omega_{S}\right)\right)^{2}\right)\left(\left(\frac{\gamma_{S}}{2}\right)^{2}+\left(\omega-\left(\omega_{S}+m_{S} \Delta \omega_{S}\right)\right)^{2}\right)} \\ & \times \frac{\left(\frac{\gamma_{I}}{2}\right)^{2}+\left(\omega-\left(\omega_{S}-m_{I}^{\prime} \Delta \omega_{I}\right)\right)\left(\omega-\left(\omega_{S}-m_{I} \Delta \omega_{I}\right)\right)}{\left(\left(\frac{\gamma_{I}}{2}\right)^{2}+\left(\omega-\left(\omega_{S}-m_{I}^{\prime} \Delta \omega_{I}\right)\right)^{2}\right)\left(\left(\frac{\gamma_{I}}{2}\right)^{2}+\left(\omega-\left(\omega_{S}-m_{I} \Delta \omega_{I}\right)\right)^{2}\right)} \end{aligned} \end{aligned}

In this paper it can be seen plotted:

Here's my first attempt to plot this, basically just by directly summing these terms (for some fairly small "truncated" amount of {ms, mi, ms', mi'}):

GHz = 10^9;
MHz = 10^6;
ps = 10^-12;
ns = 10^-9;
rulesforPlot1 = {γs -> γi, Δωi -> Δωs, τ0 -> 133 ps, Δωs -> 2 π*1.5 GHz,
γi -> 2 π*3 MHz, ωs -> 0 };

peaks = 10;
spectrumtoBeSummed =
Sinc[τ0/2 (ω - ωs)]^2 ((γs/2)^2 + (ω - (ωs + ms1 Δωs)) (ω - (ωs + ms Δωs)))/
(((γs/2)^2 + (ω - (ωs + ms1 Δωs))^2) ((γs/2)^2 +
(ω - (ωs + ms Δωs))^2)) ((γi/2)^2 +
(ω - (ωs - mi1 Δωi)) (ω - (ωs - mi Δωi)))/(((γi/2)^2 +
(ω - (ωs - ms1 Δωi))^2) ((γi/2)^2 + (ω - (ωs - mi Δωi))^2)) //. rulesforPlot1 ;
fω = Sum[spectrumtoBeSummed, {ms, -peaks, peaks}, {mi, -peaks, peaks},
{ms1, -peaks, peaks}, {mi1, -peaks, peaks}];
Plot[fω^2, {ω, -100 GHz, 100 GHz}, PlotRange -> All]


This hangs even for a small sum. I thought it would help to try to make it "numeric" instead of functional, so I tried a vectorized alternative, and it's not much better. Here is the code:

start = 0;
finish = 100 GHz;
points = 20;

peaks = 100
τvector =  N@Subdivide[start, finish, points];
ans = Flatten[spectrumtoBeSummed /. ω -> # & /@ {τvector}];
ans2 = NSum[#, {ms, -peaks, peaks}, {mi, -peaks, peaks},
{ms1, -peaks, peaks}, {mi1, -peaks, peaks} ] & /@ ans;
ListPlot[Transpose[{τvector, Flatten@ans2^2}], PlotRange -> All]


Any ideas what can be done to speed this up?

• Sorry, I don't see any modulus in your code.Do you mean the imaginary unit by $i$? Oct 30, 2020 at 7:59
• Where are $h$, $\eta$, $\gamma_1$, etc.? Oct 30, 2020 at 10:55
• @user64494, sorry, the post wasn't as clear as I thought. I made things much more clear. (In my code I used a simplified verison of the expression I linked.) I wrote out the function exactly that I want to plot, and also included an example plot in the paper I am using as a reference. Oct 30, 2020 at 21:40
• @anderstood, I updated the plots to make the constants (and the function) more clear. Sorry for it being unclear in the first place. Oct 30, 2020 at 21:40
• Also, I can't find the stackexchange link that shows how to convert my code to be more readable on this website. Oct 30, 2020 at 21:49

p0 = Plot[Sinc[x/10]^2, {x, -100, 100}, PlotStyle -> {Red, Dashed}]
p1 = ListPlot[Table[{k, Sinc[k/10]^2}, {k, -100, 100, 1.25}],
Filling -> Axis, FillingStyle -> {Darker@Gray},
Frame -> True]
Show[p0, p1]


Why is this so easy?

The author of Your source is an ambitious man and presented to You the infinite sum of profiles. That has some implicit comforts. For example, it is normed by design. This is a limit and the border has a description too. So why make the calculation more complicated than necessary.

My solution is not lacy. I assume the line profiles are small compared to the line with the envelope. The individual's line profiles are lines and the envelop is the Sinc^2.

Then I optimized to calculations time and representation comfort. I think the result is really close to the picture of the question.

Show[p0, p1, Frame -> True, FrameTicks -> {-100 + 25*Range[8], Automatic},
FrameStyle -> Thick, FrameLabel ->
{Text[Style["Signal Frequency ω-\!$$\*SubscriptBox[\(ω$$, $$S$$]\) (GHz)", 16]],
Text[Style["Normalized Signal Spectrum (a.u.)", 16]]}]


Some underpinning math:

GHz = 10^9;
MHz = 10^6;
ps = 10^-12;
ns = 10^-9;


Mhz/Ghz = 1000 = 10^3 is large enough.

NumberLinePlot[{1, 1000}]


tau = 133 ps
Δωs -> 2 π*1.5 GHz


I draw advantage from the fact that the Sinc suffices the requirement: the value of 0 is 1, the breadth is that of the dimensionless arguments always.

• Mhz/Ghz = 1000? typo? Oct 31, 2020 at 11:21