minor modif on script

maoppy/example/numerical_vs_theoretical_sigma2.py

@@ -28,9 +28,9 @@ amp = 1.0 # Phase PSD Moffat amplitude [rad²]
 alpha = np.logspace(-3,0,num=50) # Phase PSD Moffat alpha [1/m]
 ratio = 1.0
 theta = 0.0
-beta = 1.6 # Phase PSD Moffat beta power law
+beta = [1.4,1.5,1.6] # Phase PSD Moffat beta power law
 
-sig2num = np.zeros_like(alpha)
+sig2num = np.zeros((len(alpha),len(beta)))
 
 f2D = Pmodel._fxfy
 F2 = f2D[0]**2. + f2D[1]**2.
@@ -39,16 +39,19 @@ df = Pmodel._pix2freq
 mask = F2 < Fao**2
 
 for i in range(len(alpha)):
-    param = [r0,b,amp,alpha[i],ratio,theta,beta]
+    for j in range(len(beta)):
+        param = [r0,b,amp,alpha[i],ratio,theta,beta[j]]
         psd,integral = Pmodel.psd(param)
-    sig2num[i] = np.sum(psd*mask)*df**2
+        sig2num[i,j] = np.sum(psd*mask)*df**2
 
 #%% PLOT
 plt.figure(1)
 plt.clf()
-plt.semilogx(alpha/df,sig2num,label='numerical sum')
+for j in range(len(beta)):
+    plt.semilogx(alpha/df,sig2num[:,j],label='numerical sum ($\\beta=$%.1f)'%beta[j])
 plt.semilogx(alpha/df,alpha*0 + amp, label='required')
 plt.grid(which='both')
 plt.xlabel("$\\alpha$/df")
 plt.ylabel("$\\sigma^2$ [rad²]")
 plt.legend()
+plt.ylim(0,1.2)
\ No newline at end of file