How can I draw a graph?


I have a big problem.
I want to draw a graph (transmission vs intensity), how can I do?
I use this code to make a variable light source:
t = linspace(0,4000e-15, 4000);
Imin = 0;
Imax = 2.5;
nt = 4000;
I = linspace(Imin,Imax,nt);
amp = sqrt(sqrt(mu0/eps0)I10^10/(0.51.4524));
phase = -w_center

set the signal for the desired source

and I use two monitor, one of them is field time and another is frequency-domain field and power.
thank you so much,


Hi @moslem.teymoori,

If you want to get the transmission v.s. intensity graph, you could sweep through the intensity and get the transmission for each of the intensity points, then you can easily plot the figure. However, how to implement this still depends on how you define your “transmission” and “intensity”. I suppose by “intensity” you mean the intensity of your light source and according to the figure you attached, it seems to be a non-linear system. Then as you define the source, you could get its intensity as |E|^2. If you have a non-linear system, as you increase the source intensity, the non-linear part in the transmission may vary. We generally define the transmission already as normalized to the source power (also as indicated in your second figure), but in your case, if you want the un-normalized transmission, you could change the “Normalization state” to “No normalization” in the “Setting” tab. Please make sure the transmission calculated as so matches your definition of the transmission.

I hope this is want you want, otherwise, let me know your detailed requirement and share me the files :slight_smile:


I’m really thankful, but how can I sweep through the intensity and get the transmission for each of the intensity points. I worked on nonlinear materials and I mean the intensity of light source. how can I get its intensity as |E|^2.
I was wondering if you could write the script for me.
thank you so much,
test1 .lsf (435 Bytes)
test1.fsp (109.8 KB)


Hi @moslem.teymoori,

Based on your script, it seems like you want to plot this according to the time, then you don’t need to do the sweep in this case. You already have your source amplitude defined in the script (based on different time point), so you can get the light intensity from there. For the transmission, you can get it from the time monitor, the following scripts could be used:

P = getresult(“monitor”, “P”);
Pz = P.Pz;
Pz = pinch(Pz);
tP = P.t;
Pt = matrix(length(t));
for (ii = 1:length(t)) {
index = find(tP,t(ii));
Pt(ii) = Pz(index);

and Pz is the power component in z-dimentsion. The for loop is used to pick out the transmission data that matches your time points defined for the source.

I hope this could help :slight_smile:


Thank you so much, it’s totally useful but I work on optical bistability and I want that the plot become normalized such as this:



For optical bistability, you can normalize the results using a similar method to the bistability_results.lsf script file from this example:

The example uses a time monitor in front of the source to measure the injected fields over time, and another time monitor which measures the transmitted fields over time. A simple normalized result would be the measured transmitted field intensity divided by the injected field intensity at each point in time. However, it turns out that there is some high frequency noise in the resulting the simply normalized result, so an additional high frequency filter is applied in the script to get the final normalized result.

Since the intensity is increased linearly over time, the normalized transmission over time can be related to the normalized transmission over input intensity.

As an example, I modified your simulation file by adding a time monitor in front of the source to measure the input intensity, and I changed the polarization of the source to Ey. I then used a script file which I modified from the bistability_results.lsf script file.
modified_test1.fsp (111.9 KB)
modified_bistability_results.lsf (815 Bytes)

If you try it out you will see that the result is close to 1 since the material is not nonlinear right now, however the start and end range of the plot do have some error which is due to the high frequency filter.

Hopefully this helps!


Thank you so much,
I have a question.
If I use the polarization of 45 degrees, what should I do? Can I change Ey to power in z direction?
I’m really thankful.


If you use a 45 degree polarization angle, the total E field will have both Ex and Ey components. You can get the total E field amplitude by collecting both Ex and Ey from the monitor, then E_total = sqrt(abs(Ex)^2+abs(Ey)^2).

In the script file, you would need to replace Ey with E_total to get the total normalized transmission. If you need further help with this, let me know!


sorry, I didn’t understand.
Is this code true:
w0 = 2pic/1482e-9; # center frequency
Imin = 10; Imax = 10;

Ey = getdata(‘monitor’,‘Ey’);
Esy = getdata(‘monitor_7’,‘Ey’);
t = getdata(‘monitor’,‘t’);
I = linspace(Imin, Imax, length(t));
zero_pad = 2^16;
w = fftw(t,1,zero_pad);
filter1 = 2*exp(-(w-w0)^2/(150e12)^2);

Ex = getdata(‘monitor’,‘Ex’);
Esx = getdata(‘monitor_7’,‘Ex’);
t = getdata(‘monitor’,‘t’);
I = linspace(Imin, Imax, length(t));
zero_pad = 2^16;
w = fftw(t,1,zero_pad);
filter1 = 2*exp(-(w-w0)^2/(150e12)^2);

remove high frequency components y

Ey_through_w = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Ey),1,zero_pad);
Ey_through_t = invfft(pinch(Ey_through_w)filter1);
Ey_through_t = Ey_through_t(1:length(t));
Esy_through_w = 2
Esy_through_t = invfft(pinch(Esy_through_w)*filter1);
Esy_through_t = Esy_through_t(1:length(t));

remove high frequency components x

Ex_through_w = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Ey),1,zero_pad);
Ex_through_t = invfft(pinch(Ey_through_w)filter1);
Ex_through_t = Ey_through_t(1:length(t));
Esx_through_w = 2
Esx_through_t = invfft(pinch(Esy_through_w)*filter1);
Esx_through_t = Esy_through_t(1:length(t));

plot results

plot(I, sqrt(abs(Ex)^2+abs(Ey)^2)/sqrt(abs(Esx)^2+abs(Esy)^2),‘I (MW/cm^2)’,’(E_out/E_source)^2’);
legend(“Increasing Source Intensity”);


Hi @moslem.teymoori,

I think what @nlui means is that yo now have to consider the x-component of the field. So you need to read out the x-component field from the monitor too. Basically you just need to repeat the same thing as you did to the y component. e.g., add the following scripts to the original modified_bistability_results.lsf file:

Ex = getdata(‘monitor’,‘Ex’);
Esx = getdata(‘monitor_1’,‘Ex’);

Ex_through_w = 2*((1:length(w))<=(length(w)/2+0.1))*fft(pinch(Ex),1,zero_pad);
Ex_through_t = invfft(pinch(Ex_through_w)filter1);
Ex_through_t = Ex_through_t(1:length(t));
Esx_through_w = 2
Esx_through_t = invfft(pinch(Esx_through_w)*filter1);
Esx_through_t = Esx_through_t(1:length(t));

E_total = sqrt(abs(Ex_through_t )^2+abs(Ey_through_t )^2);
Es_total = sqrt(abs(Esx_through_t )^2+abs(Esy_through_t )^2);

plot(I, E_total^2/Es_total^2, ‘I (MW/cm^2)’ , ‘(E_out/E_source)^2’);
legend(“Increasing Source Intensity”);

Hope this helps :slight_smile:


thank you so much,
it’s really useful.
I have a question, what is the role of “remove high frequency components y or x” exactly?


Hi @moslem.teymoori,

The fft command with option1 (which is the second argument of the function, please refer to the online help page: as 1 performs the the standard fft with zero frequency at the first element of the matrix. Then the “remove high frequency components” part removes the high frequency component from the matrix and leaves only the result up to the Nyquist frequency.

In your case, since your time signal is all real, you could also use the fft command option1 as 2, which zero frequency is the first element, but only data up to and including the Nyquist frequency is stored.

I hope this could help :slight_smile: