Update mopc.py

samodeo · Nov 25, 2022 · 66cd2ff · 66cd2ff
1 parent 050d0e2
commit 66cd2ff
Showing 1 changed file with 50 additions and 44 deletions.
diff --git a/mopc/mopc.py b/mopc/mopc.py
@@ -53,7 +53,8 @@ def project_prof_beam(tht,M,z,theta,nu,f_beam):
     disc_fac = np.sqrt(2)
     l0 = 30000.
     NNR = 100
-    NNR2 = 3*NNR
+    resolution_factor = 3.0
+    NNR2 = resolution_factor*NNR
 
     #fwhm = beam
     #fwhm *= np.pi / (180.*60.)
@@ -88,8 +89,8 @@ def project_prof_beam(tht,M,z,theta,nu,f_beam):
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta_smooth = thta_smooth[:,None]
     thta2_smooth = thta2_smooth[:,None]
@@ -102,9 +103,9 @@ def project_prof_beam(tht,M,z,theta,nu,f_beam):
     Pth2D = 2*np.trapz(Pth(rint,M,z,theta_p), x=rad * kpc_cgs, axis=1) * 1e3
     Pth2D2 = 2*np.trapz(Pth(rint2,M,z,theta_p), x=rad2 * kpc_cgs, axis=1) * 1e3
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta = thta[:,None,None]
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2 = thta2[:,None,None]
 
     phi = np.linspace(0., 2*np.pi, 100)
@@ -130,8 +131,8 @@ def project_prof_beam(tht,M,z,theta,nu,f_beam):
             * f_beam(np.sqrt(thta2**2+thta2_smooth**2-2*thta2*thta2_smooth*np.cos(phi))) , x=phi,axis=2)
 
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     rho2D_beam = np.trapz(rho2D_beam0, x=thta_smooth,axis=1)
     rho2D2_beam = np.trapz(rho2D2_beam0, x=thta2_smooth,axis=1)
@@ -160,7 +161,8 @@ def project_prof_beam_y(tht,M,z,theta,f_beam,y_beam):
     disc_fac = np.sqrt(2)
     l0 = 30000.
     NNR = 100
-    NNR2 = 2*NNR
+    resolution_factor = 2.0
+    NNR2 = resolution_factor*NNR
 
     fwhm = y_beam
     fwhm *= np.pi / (180.*60.)
@@ -191,8 +193,8 @@ def project_prof_beam_y(tht,M,z,theta,f_beam,y_beam):
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta_smooth = thta_smooth[:,None]
     thta2_smooth = thta2_smooth[:,None]
@@ -207,9 +209,9 @@ def project_prof_beam_y(tht,M,z,theta,f_beam,y_beam):
     Pth2D2 = 2*np.trapz(Pth(rint2,M,z,theta_p), x=rad2 * kpc_cgs, axis=1) * 1e3
 
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta = thta[:,None,None]
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2 = thta2[:,None,None]
 
     phi = np.linspace(0., 2*np.pi, 100)
@@ -222,9 +224,9 @@ def project_prof_beam_y(tht,M,z,theta,f_beam,y_beam):
 
     thta_y = (np.arange(NNR) + 1.)*dtht
     thta2_y = (np.arange(NNR) + 1.)*dtht2
-    thta_smooth_y = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth_y = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta_y = thta_y[:,None]
-    thta2_smooth_y = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth_y = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2_y = thta2_y[:,None]
 
     rho2D_beam0 = np.trapz( thta_smooth  * rho2D  
@@ -239,8 +241,8 @@ def project_prof_beam_y(tht,M,z,theta,f_beam,y_beam):
                                  * special.iv(0, thta2_smooth_y*thta2_y/ sigmaBeam**2), x=thta2_smooth_y, axis=1)
 
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
@@ -271,7 +273,8 @@ def project_prof_beam_sim_rho(tht,M,z,theta_rho,f_beam):
     disc_fac = np.sqrt(2)
     l0 = 30000.
     NNR = 100
-    NNR2 = 3.*NNR
+    resolution_factor = 3.0
+    NNR2 = resolution_factor*NNR
 
     #fwhm = beam #arcmin
     #fwhm *= np.pi / (180.*60.) #rad
@@ -297,8 +300,8 @@ def project_prof_beam_sim_rho(tht,M,z,theta_rho,f_beam):
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta_smooth = thta_smooth[:,None]
     thta2_smooth = thta2_smooth[:,None]
@@ -309,9 +312,9 @@ def project_prof_beam_sim_rho(tht,M,z,theta_rho,f_beam):
     rho2D = 2*np.trapz(rho_gnfw(rint,M,z,theta_sim_rho), x=rad * kpc_cgs, axis=1) * 1e3
     rho2D2 = 2*np.trapz(rho_gnfw(rint2,M,z,theta_sim_rho), x=rad2 * kpc_cgs, axis=1) * 1e3
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta = thta[:,None,None]
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2 = thta2[:,None,None]
 
     phi = np.linspace(0., 2*np.pi, 100)
@@ -328,8 +331,8 @@ def project_prof_beam_sim_rho(tht,M,z,theta_rho,f_beam):
     rho2D2_beam0 = np.trapz( thta2_smooth  * rho2D2  
             * f_beam(np.sqrt(thta2**2+thta2_smooth**2-2*thta2*thta2_smooth*np.cos(phi))) , x=phi,axis=2)
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     rho2D_beam = np.trapz(rho2D_beam0, x=thta_smooth,axis=1)
     rho2D2_beam = np.trapz(rho2D2_beam0, x=thta2_smooth,axis=1)
@@ -342,7 +345,7 @@ def project_prof_beam_sim_rho(tht,M,z,theta_rho,f_beam):
     sig  = 2.0*np.pi*dtht *np.sum(thta *rho2D_beam)
     sig2 = 2.0*np.pi*dtht2*np.sum(thta2*rho2D2_beam)
 
-    sig_all_beam = (2*sig - sig2) * v_rms * ST_CGS * TCMB * 1e6 * ((2. + 2.*XH)/(3.+5.*XH)) / MP_CGS #/ (np.pi * np.radians(tht/60.)**2)
+    sig_all_beam = (2*sig - sig2) * v_rms * ST_CGS * TCMB * 1e6 * ((1. + XH)/2.) / MP_CGS #/ (np.pi * np.radians(tht/60.)**2)
 
     return sig_all_beam
 
@@ -351,7 +354,8 @@ def project_prof_beam_rho_dm(tht,M,z,f_beam):
     disc_fac = np.sqrt(2)
     l0 = 30000.
     NNR = 100
-    NNR2 = 3.*NNR
+    resolution_factor = 3.0
+    NNR2 = resolution_factor*NNR
 
     #fwhm = beam #arcmin
     #fwhm *= np.pi / (180.*60.) #rad
@@ -380,8 +384,8 @@ def project_prof_beam_rho_dm(tht,M,z,f_beam):
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta_smooth = thta_smooth[:,None]
     thta2_smooth = thta2_smooth[:,None]
@@ -392,9 +396,9 @@ def project_prof_beam_rho_dm(tht,M,z,f_beam):
     rho2D = 2*np.trapz(rho_dm(rint,M,z), x=rad * kpc_cgs, axis=1) * 1e3
     rho2D2 = 2*np.trapz(rho_dm(rint2,M,z), x=rad2 * kpc_cgs, axis=1) * 1e3
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta = thta[:,None,None]
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2 = thta2[:,None,None]
 
     phi = np.linspace(0., 2*np.pi, 100)
@@ -411,8 +415,8 @@ def project_prof_beam_rho_dm(tht,M,z,f_beam):
     rho2D2_beam0 = np.trapz( thta2_smooth  * rho2D2  
             * f_beam(np.sqrt(thta2**2+thta2_smooth**2-2*thta2*thta2_smooth*np.cos(phi))) , x=phi,axis=2)
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     rho2D_beam = np.trapz(rho2D_beam0, x=thta_smooth,axis=1)
     rho2D2_beam = np.trapz(rho2D2_beam0, x=thta2_smooth,axis=1)
@@ -426,7 +430,7 @@ def project_prof_beam_rho_dm(tht,M,z,f_beam):
     sig  = 2.0*np.pi*dtht *np.sum(thta *rho2D_beam)
     sig2 = 2.0*np.pi*dtht2*np.sum(thta2*rho2D2_beam)
 
-    sig_all_beam = (2*sig - sig2) * v_rms * ST_CGS * TCMB * 1e6 * ((2. + 2.*XH)/(3.+5.*XH)) / MP_CGS #/ (np.pi * np.radians(tht/60.)**2)
+    sig_all_beam = (2*sig - sig2) * v_rms * ST_CGS * TCMB * 1e6 * ((1. + XH)/2.) / MP_CGS #/ (np.pi * np.radians(tht/60.)**2)
 
     return sig_all_beam
 
@@ -436,7 +440,8 @@ def project_prof_beam_sim_pth(tht,M,z,theta_pth,nu,f_beam):
     disc_fac = np.sqrt(2)
     l0 = 30000.
     NNR = 100
-    NNR2 = 3.5*NNR
+    resolution_factor = 3.5
+    NNR2 = resolution_factor*NNR
 
     #fwhm = beam
     #fwhm *= np.pi / (180.*60.)
@@ -468,8 +473,8 @@ def project_prof_beam_sim_pth(tht,M,z,theta_pth,nu,f_beam):
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta_smooth = thta_smooth[:,None]
     thta2_smooth = thta2_smooth[:,None]
@@ -480,9 +485,9 @@ def project_prof_beam_sim_pth(tht,M,z,theta_pth,nu,f_beam):
     Pth2D = 2*np.trapz(Pth_gnfw(rint,M,z,theta_sim_pth), x=rad * kpc_cgs, axis=1) * 1e3
     Pth2D2 = 2*np.trapz(Pth_gnfw(rint2,M,z,theta_sim_pth), x=rad2 * kpc_cgs, axis=1) * 1e3
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta = thta[:,None,None]
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2 = thta2[:,None,None]
     phi = np.linspace(0., 2*np.pi, 50)
     phi = phi[None,None,:]
@@ -497,8 +502,8 @@ def project_prof_beam_sim_pth(tht,M,z,theta_pth,nu,f_beam):
     Pth2D2_beam0 = np.trapz( thta2_smooth  * Pth2D2  
             * f_beam(np.sqrt(thta2**2+thta2_smooth**2-2*thta2*thta2_smooth*np.cos(phi))) , x=phi,axis=2)
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     Pth2D_beam = np.trapz(Pth2D_beam0, x=thta_smooth,axis=1)
     Pth2D2_beam = np.trapz(Pth2D2_beam0, x=thta2_smooth,axis=1)
@@ -522,7 +527,8 @@ def project_prof_beam_sim_y(tht,M,z,theta_pth, beam):
     disc_fac = np.sqrt(2)
     l0 = 30000.
     NNR = 100
-    NNR2 = 3.5*NNR
+    resolution_factor = 3.5
+    NNR2 = resolution_factor*NNR
 
     fwhm = beam
     fwhm *= np.pi / (180.*60.)
@@ -552,8 +558,8 @@ def project_prof_beam_sim_y(tht,M,z,theta_pth, beam):
     thta = (np.arange(NNR) + 1.)*dtht
     thta2 = (np.arange(NNR) + 1.)*dtht2
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
 
     thta_smooth = thta_smooth[:,None]
     thta2_smooth = thta2_smooth[:,None]
@@ -564,9 +570,9 @@ def project_prof_beam_sim_y(tht,M,z,theta_pth, beam):
     Pth2D = 2*np.trapz(Pth_gnfw(rint,M,z,theta_sim_pth), x=rad * kpc_cgs, axis=1) * 1e3
     Pth2D2 = 2*np.trapz(Pth_gnfw(rint2,M,z,theta_sim_pth), x=rad2 * kpc_cgs, axis=1) * 1e3
 
-    thta_smooth = (np.arange(NNR2) + 1.)*dtht
+    thta_smooth = (np.arange(NNR2) + 1.)*dtht/resolution_factor
     thta = thta[:,None]
-    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2
+    thta2_smooth = (np.arange(NNR2) + 1.)*dtht2/resolution_factor
     thta2 = thta2[:,None]
 
     Pth2D_beam  = np.trapz(thta_smooth  * Pth2D  * np.exp(-0.5*thta_smooth**2 /sigmaBeam**2)