Python Planck13.comoving_distance方法代碼示例

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def slice_z_cosmo(slice, z, radius):
    """Slices in redshift considering a distance to cut"""
    zmax = z_at_value(cosmo.comoving_distance, cosmo.comoving_distance(z) + radius)
    zmin = z_at_value(cosmo.comoving_distance, cosmo.comoving_distance(z) - radius)

    print zmin - zmax

    slice = slice[np.where(slice["zb_1"] < zmax)]
    slice = slice[np.where(slice["zb_1"] > zmin)]
    return slice

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def single_run_test(ind,ysc1,ysc2,q,vd,pha,zl,zs):
    dsx_sdss     = 0.396         # pixel size of SDSS detector.
    R  = 3.0000     #
    nnn = 300      #Image dimension
    bsz = 9.0 # arcsecs
    dsx = bsz/nnn         # pixel size of SDSS detector.
    nstd = 59 #^2

    xx01 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xx02 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xi2,xi1 = np.meshgrid(xx01,xx02)
    #----------------------------------------------------------------------
    dsi = 0.03
    g_source = pyfits.getdata("./439.0_149.482739_1.889989_processed.fits")
    g_source = np.array(g_source,dtype="<d")*10.0
    g_source[g_source<=0.0001] = 1e-6
    #----------------------------------------------------------------------
    xc1 = 0.0       #x coordinate of the center of lens (in units of Einstein radius).
    xc2 = 0.0       #y coordinate of the center of lens (in units of Einstein radius).
    #q   = 0.7       #Ellipticity of lens.
    rc  = 0.0       #Core size of lens (in units of Einstein radius).
    re  = re_sv(vd,zl,zs)       #Einstein radius of lens.
    #pha = 45.0      #Orintation of lens.
    lpar = np.asarray([xc1,xc2,q,rc,re,pha])
    #----------------------------------------------------------------------
    ai1,ai2,mua = lens_equation_sie(xi1,xi2,lpar)

    yi1 = xi1-ai1
    yi2 = xi2-ai2

    g_limage = lv4.call_ray_tracing(g_source,xi1,xi2,ysc1,ysc2,dsi)
    g_limage[g_limage<=0.0001] = 1e-6
    g_limage = p2p.cosccd2mag(g_limage)
    g_limage = p2p.mag2sdssccd(g_limage)

    #-------------------------------------------------------------
    # Need to be Caliborate the mags
    dA = Planck13.comoving_distance(zl).value*1000./(1+zl)
    Re = dA*np.sin(R*np.pi/180./3600.)
    counts  =Brightness(Re,vd)
    vpar = np.asarray([counts,R,xc1,xc2,q,pha])
    #g_lens = deVaucouleurs(xi1,xi2,xc1,xc2,counts,R,1.0-q,pha)
    g_lens = de_vaucouleurs_2d(xi1,xi2,vpar)

    #pl.figure()
    #pl.contourf(g_lens)
    #pl.colorbar()

    g_clean_ccd1 = g_lens*0.0+g_limage
    g_clean_ccd2 = g_lens*1.0+g_limage
    from scipy.ndimage.filters import gaussian_filter
    pl.figure()
    pl.contourf(g_clean_ccd1)
    pl.colorbar()
    pl.savefig("/Users/uranus/GitHub/data_arc_finding_cnn/unlensed_output_pngs/"+str(i)+"_lensed_imgs_only.png")
    #-------------------------------------------------------------
    g_images_psf1 = gaussian_filter(g_clean_ccd1, 2.0)
    g_images_psf2 = gaussian_filter(g_clean_ccd2, 2.0)
    #g_images_psf = ss.convolve(g_clean_ccd,g_psf,mode="same")
    #g_images_psf = g_clean_ccd
    #-------------------------------------------------------------
    # Need to be Caliborate the mags
    g_noise = noise_map(nnn,nnn,np.sqrt(nstd),"Gaussian")
    #output_filename = "./output_fits/noise_map.fits"
    #pyfits.writeto(output_filename,g_noise,clobber=True)
    #-------------------------------------------------------------
    g_final1 = g_images_psf1+g_noise
    output_filename = "/Users/uranus/GitHub/data_arc_finding_cnn/unlensed_output_fits/"+str(i)+"_lensed_imgs_only.fits"
    pyfits.writeto(output_filename,g_final1,clobber=True)
    #-------------------------------------------------------------
    g_final2 = g_images_psf2+g_noise
    output_filename = "/Users/uranus/GitHub/data_arc_finding_cnn/unlensed_output_fits/"+str(i)+"_all_imgs.fits"
    pyfits.writeto(output_filename,g_final2,clobber=True)

    return 0

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def single_run_test(ind, ysc1, ysc2, q, vd, pha, zl, zs, lens_tag=1):
    # dsx_sdss     = 0.396         # pixel size of SDSS detector.
    R = 3.0
    nnn = 400  # Image dimension
    bsz = 9.0  # arcsecs
    dsx = bsz / nnn         # pixel size of SDSS detector.
    nstd = 59  # ^2

    xi1, xi2 = make_r_coor(nnn, dsx)
# ----------------------------------------------------------------------
    dsi = 0.03
    g_source = pyfits.getdata(
        "./gals_sources/439.0_149.482739_1.889989_processed.fits")
    g_source = np.array(g_source, dtype="<d") * 10.0
    g_source[g_source <= 0.0001] = 1e-6
# ----------------------------------------------------------------------
    # x coordinate of the center of lens (in units of Einstein radius).
    xc1 = 0.0
    # y coordinate of the center of lens (in units of Einstein radius).
    xc2 = 0.0
    rc = 0.0  # Core size of lens (in units of Einstein radius).
    re = re_sv(vd, zl, zs)  # Einstein radius of lens.
    print "re = ", re
    re_sub = 0.0 * re
    a_sub = a_b_bh(re_sub, re)
    ext_shears = 0.06
    ext_angle = -0.39
    ext_kappa = 0.08

    # ext_shears = 0.0
    # ext_angle = 0.0
    # ext_kappa = 0.0

# ----------------------------------------------------------------------
    #lpar = np.asarray([xc1, xc2, q, rc, re, pha])
    #ai1, ai2, kappa_out, shr1, shr2, mua = lensing_signals_sie(xi1, xi2, lpar)
    #ar = np.sqrt(ai1 * ai1 + ai2 * ai2)

    # psi_nie = psk.potential_nie(xc1, xc2, pha, q, re, rc, ext_shears, ext_angle,
    #                            ext_kappa, xi1, xi2)
    #ai1, ai2 = np.gradient(psi_nie, dsx)

    ai1, ai2 = psk.deflection_nie(xc1, xc2, pha, q, re, rc, ext_shears, ext_angle,
                                  ext_kappa, xi1, xi2)

    as1, as2 = psk.deflection_sub_pJaffe(0.0, -2.169, re_sub, 0.0, a_sub, xi1, xi2)

    yi1 = xi1 - ai1 - as1
    yi2 = xi2 - ai2 - as2

    g_limage = lv4.call_ray_tracing(g_source, yi1, yi2, ysc1, ysc2, dsi)
    g_limage[g_limage <= 0.25] = 1e-6

    # pl.figure()
    # pl.contourf(g_limage)
    # pl.colorbar()

    g_limage = cosccd2mag(g_limage)
    g_limage = mag2sdssccd(g_limage)

    # pl.figure()
    # pl.contourf(g_limage)
    # pl.colorbar()
# -------------------------------------------------------------
    dA = Planck13.comoving_distance(zl).value * 1000. / (1.0 + zl)
    Re = dA * np.sin(R * np.pi / 180. / 3600.)
    counts = Brightness(Re, vd)
    vpar = np.asarray([counts, R, xc1, xc2, q, pha])
    g_lens = de_vaucouleurs_2d(xi1, xi2, vpar)

    g_clean_ccd = g_lens * lens_tag + g_limage
    output_filename = "./fits_outputs/clean_lensed_imgs.fits"
    pyfits.writeto(output_filename, g_clean_ccd, overwrite=True)
# -------------------------------------------------------------
    from scipy.ndimage.filters import gaussian_filter
    g_images_psf = gaussian_filter(g_clean_ccd, 2.0)
# -------------------------------------------------------------
    g_noise = noise_map(nnn, nnn, np.sqrt(nstd), "Gaussian")
    output_filename = "./fits_outputs/noise_map.fits"
    pyfits.writeto(output_filename, g_noise, overwrite=True)
    g_final = g_images_psf + g_noise
# -------------------------------------------------------------
    g_clean_ccd = g_limage
    g_images_psf = gaussian_filter(g_clean_ccd, 2.0)
    g_final = g_images_psf + g_noise
    output_filename = "./fits_outputs/lensed_imgs_only.fits"
    pyfits.writeto(output_filename, g_final, overwrite=True)
# -------------------------------------------------------------
    output_filename = "./fits_outputs/full_lensed_imgs.fits"
    pyfits.writeto(output_filename, g_final + g_lens, overwrite=True)

    # pl.figure()
    # pl.contourf(g_final)
    # pl.colorbar()

    return 0

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
#Find the RA, DEC, Z
#galra=galdat.RA[:100]
#galdec=galdat.DEC[:100]
#galz=galdat.Z[:100]

cutra=cutdat.RA
cutdec=cutdat.DEC
cutz=cutdat.Z
randra=randat.RA
randdec=randat.DEC
randz=randat.Z


#Convert redshift to Mpc (Using Planck Results stored in astropy.cosmology)
#comoving=Planck13.comoving_distance(galz)
cutX=Planck13.comoving_distance(cutz)
randX=Planck13.comoving_distance(randz)
end2=time.time()
print end2-start, 'Comoving distances calculated'

'''
#From working cosmology code
def Horizon(x, t, Om, Og, Ol, Ok):
	F=(x**-2)*(Om/x**3+Og/x**4+Ol+Ok/x**2)**-0.5
	return F
Om=0.307
Og=0
Ol=0.693
Ok=1-(Om+Og+Ol)
t=linspace(0,1,100000)
DA=[]

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def single_run_test(ind,ysc1,ysc2,q,vd,pha,zl,zs):
    dsx_sdss     = 0.396         # pixel size of SDSS detector.
    R  = 3.0000     #
    #zl = 0.2     #zl is the redshift of the lens galaxy.
    #zs = 1.0
    #vd = 520    #Velocity Dispersion.
    nnn = 128      #Image dimension
    bsz = dsx_sdss*nnn # arcsecs
    dsx = bsz/nnn         # pixel size of SDSS detector.
    nstd = 59 #^2

    xx01 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xx02 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xi2,xi1 = np.meshgrid(xx01,xx02)
    #----------------------------------------------------------------------
    #ysc1 = 0.2
    #ysc2 = 0.5
    dsi = 0.03
    g_source = pyfits.getdata("./439.0_149.482739_1.889989_processed.fits")
    g_source = np.array(g_source,dtype="<d")*10.0
    g_source[g_source<=0.0001] = 1e-6
    #print np.sum(g_source)
    #print np.max(g_source)
    #pl.figure()
    #pl.contourf(g_source)
    #pl.colorbar()
    #g_source = p2p.cosccd2mag(g_source)
    ##g_source = p2p.mag2sdssccd(g_source)
    ##print np.max(g_source*13*13*52.0)
    #pl.figure()
    #pl.contourf(g_source)
    #pl.colorbar()
    #----------------------------------------------------------------------
    xc1 = 0.0       #x coordinate of the center of lens (in units of Einstein radius).
    xc2 = 0.0       #y coordinate of the center of lens (in units of Einstein radius).
    #q   = 0.7       #Ellipticity of lens.
    rc  = 0.0       #Core size of lens (in units of Einstein radius).
    re  = re_sv(vd,zl,zs)       #Einstein radius of lens.
    #pha = 45.0      #Orintation of lens.
    lpar = np.asarray([xc1,xc2,q,rc,re,pha])
    #----------------------------------------------------------------------
    ai1,ai2,mua = lens_equation_sie(xi1,xi2,lpar)

    yi1 = xi1-ai1
    yi2 = xi2-ai2

    g_limage = lv4.call_ray_tracing(g_source,yi1,yi2,ysc1,ysc2,dsi)
    g_limage[g_limage<=0.0001] = 1e-6
    g_limage = p2p.cosccd2mag(g_limage)
    g_limage = p2p.mag2sdssccd(g_limage)

    #pl.figure()
    #pl.imshow((g_limage),interpolation='lanczos',cmap=cm.gray)
    #pl.colorbar()

    #-------------------------------------------------------------
    # Need to be Caliborate the mags
    dA = Planck13.comoving_distance(zl).value*1000./(1+zl)
    Re = dA*np.sin(R*np.pi/180./3600.)
    counts  =Brightness(Re,vd)
    vpar = np.asarray([counts,R,xc1,xc2,q,pha])
    #g_lens = deVaucouleurs(xi1,xi2,xc1,xc2,counts,R,1.0-q,pha)
    g_lens = de_vaucouleurs_2d(xi1,xi2,vpar)

    #pl.figure()
    #pl.imshow((g_lens),interpolation='nearest',cmap=cm.gray)
    #pl.colorbar()

    g_clean_ccd = g_lens+g_limage

    #pl.figure()
    #pl.imshow((g_clean_ccd),interpolation='nearest',cmap=cm.gray)
    #pl.colorbar()

    g_clean_ccd = congrid.congrid(g_clean_ccd,[128,128])

    #-------------------------------------------------------------
    file_psf = "../PSF_and_noise/sdsspsf.fits"
    g_psf = pyfits.getdata(file_psf)-1000.0
    g_psf = g_psf/np.sum(g_psf)

    #new_shape=[0,0]
    #new_shape[0]=np.shape(g_psf)[0]*dsx_sdss/dsx
    #new_shape[1]=np.shape(g_psf)[1]*dsx_sdss/dsx
    #g_psf = rebin_psf(g_psf,new_shape)

    g_images_psf = ss.fftconvolve(g_clean_ccd,g_psf,mode="same")
    #g_images_psf = ss.convolve(g_clean_ccd,g_psf,mode="same")
    #g_images_psf = g_clean_ccd

    #pl.figure()
    #pl.imshow((g_psf),interpolation='nearest',cmap=cm.gray)
    #pl.colorbar()

    #-------------------------------------------------------------
    # Need to be Caliborate the mags
    #g_noise = noise_map(nnn,nnn,np.sqrt(nstd),"Gaussian")
    g_noise = noise_map(128,128,np.sqrt(nstd),"Gaussian")
    g_final = g_images_psf+g_noise

#.........這裏部分代碼省略.........

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def X(f):
    z=f21/f-1
    return pl13.comoving_distance(z).value

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def kperp2u(kperp,z):
    return kperp/(2.*pi)*pl13.comoving_distance(z).value

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def u2kperp(u,z):
    return u*2.*pi/pl13.comoving_distance(z).value

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
                    help="Noise variance",metavar="value",
                    type="float")
parser.add_option("--NoiseType",dest="NoiseType",default="Gaussian",
                    help="Gaussian or Poisson",metavar="value",
                    type="string")

(o,args)=parser.parse_args()
R       =o.R_eff     #vd is velocity dispersion.
zl      =o.zlens     #zl is the redshift of the lens galaxy.
Vd      =o.VelDis    #Velocity Dispersion.
Npix    =o.Npix      #Image dimension
e       =o.E         #Ellipticity of simulated galaxy
phi     =o.Phi       #Orientation angle wrt North
nstd    =np.sqrt(o.NoiseVar)

dA      =Planck13.comoving_distance(zl).value*1000./(1+zl)
Re      =dA*np.sin(R*np.pi/180./3600.)
print Re
imgal   =deVaucouleurs(R,Re,Vd,e,phi,Npix)

skycount=sky_r/(nMgyCount_r)
if o.NoiseType=='Poisson':
        print '## Hey! You choose Poisson noise. Which is not good for now.'
        noise=np.random.poisson(nstd,(Npix,Npix))-nstd
if o.NoiseType=='Gaussian':
        print '## Hey! You choose Gaussian noise.'
        noise=nstd*np.random.normal(0.0,1.0,(Npix,Npix))


import matplotlib.pyplot as plt
#plt.imshow((imgal+noise),interpolation='Nearest',cmap=cm.gray)

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def single_run_test(ind,ysc1,ysc2,q,vd,pha,zl,zs):
    zeropoint = 18

    dsx_sdss     = 0.396         # pixel size of SDSS detector.


    R  = 2.9918     #vd is velocity dispersion.
    #zl = 0.2     #zl is the redshift of the lens galaxy.
    #zs = 1.0
    #vd = 520    #Velocity Dispersion.
    nnn = 512      #Image dimension
    bsz = 30.0 # arcsecs
    dsx = bsz/nnn         # pixel size of SDSS detector.
    nstd = 59

    xx01 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xx02 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xi2,xi1 = np.meshgrid(xx01,xx02)
    #----------------------------------------------------------------------
    #ysc1 = 0.2
    #ysc2 = 0.5
    dsi = 0.03
    g_source = pyfits.getdata("./439.0_149.482739_1.889989_processed.fits")
    g_source = np.array(g_source,dtype="<d")
    g_source = p2p.pixcos2pixsdss(g_source)
    #----------------------------------------------------------------------
    xc1 = 0.0       #x coordinate of the center of lens (in units of Einstein radius).
    xc2 = 0.0       #y coordinate of the center of lens (in units of Einstein radius).
    #q   = 0.7       #Ellipticity of lens.
    rc  = 0.0       #Core size of lens (in units of Einstein radius).
    re  = re_sv(vd,zl,zs)       #Einstein radius of lens.
    #pha = 45.0      #Orintation of lens.
    lpar = np.asarray([xc1,xc2,q,rc,re,pha])
    #----------------------------------------------------------------------
    ai1,ai2,mua = lens_equation_sie(xi1,xi2,lpar)

    yi1 = xi1-ai1
    yi2 = xi2-ai2

    g_limage = lv4.call_ray_tracing(g_source,yi1,yi2,ysc1,ysc2,dsi)
    g_limage = mag_to_flux(g_limage,zeropoint)

    #pl.figure()
    #pl.contourf(xi1,xi2,g_limage)
    #pl.colorbar()
    #-------------------------------------------------------------
    # Need to be Caliborate the mags
    dA = Planck13.comoving_distance(zl).value*1000./(1+zl)
    Re = dA*np.sin(R*np.pi/180./3600.)
    counts  =Brightness(R,vd)
    vpar = np.asarray([counts,Re,xc1,xc2,q,pha])
    #g_lens = deVaucouleurs(xi1,xi2,xc1,xc2,counts,R,1.0-q,pha)
    g_lens = de_vaucouleurs_2d(xi1,xi2,vpar)

    g_lens = ncounts_to_flux(g_lens*1.5e-4,zeropoint)
    #-------------------------------------------------------------
    file_psf = "../PSF_and_noise/sdsspsf.fits"
    g_psf = pyfits.getdata(file_psf)-1000.0
    g_psf = g_psf/np.sum(g_psf)
    new_shape=[0,0]
    new_shape[0]=np.shape(g_psf)[0]*dsx_sdss/dsx
    new_shape[1]=np.shape(g_psf)[1]*dsx_sdss/dsx
    g_psf = rebin_psf(g_psf,new_shape)
    print(np.max(g_psf))
    g_limage = ss.fftconvolve(g_limage+g_lens,g_psf,mode="same")

    #pl.figure()
    #pl.contourf(xi1,xi2,g_limage)
    #pl.colorbar()
    #-------------------------------------------------------------
    # Need to be Caliborate the mags
    g_noise = noise_map(nnn,nnn,nstd,"Gaussian")
    g_noise = ncounts_to_flux(g_noise*1e-0+skycount,zeropoint)
    g_limage = g_limage+g_noise

    print np.shape(g_limage)
    g_limage = congrid.congrid(g_limage,[128,128])
    g_limage = g_limage-np.min(g_limage)

    pl.figure()
    #pl.contourf(xi1,xi2,g_limage)
    pl.contourf(g_limage)
    pl.colorbar()
    #-------------------------------------------------------------

    output_filename = "../output_fits/"+str(ind)+".fits"
    pyfits.writeto(output_filename,g_limage,clobber=True)

    pl.show()

    return 0

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
 def distance_Mpc(self):
     return cosmo.comoving_distance(self.redshift).value

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
#---------------------------

#------ dirs --
home=os.getenv('HOME')
scratch='/global/cscratch1/sd/fantaye' #os.getenv('SCRATCH')
srcdir=join(scratch,'seedlets/fits1024/SurfDensMap/')
#os.listdir(srcdir)
print('home:',home)
print('srcdir',srcdir)


#---- original file info
#columns: ##1:snap 2:redshift 3:Nparts 4:shot_noise
id_sim, z_sim, nobj_sim,snoise_sim = np.loadtxt('../data/3Dneedlets/std_total_parts_found_from_049_to_049.dat',unpack=True)
r_sim=cosmo.comoving_distance(z_sim)
rho_file=lambda x:   '../data/3Dneedlets/HealPixMap_nside1024_NoRnd_Radius0_std/SurfDensMap_snap_049.{}.fits'.format(x)
kappa_file=lambda x: '../data/3Dneedlets/LightCones_nside1024_NoRnd_Radius0_std/KappaMap_snap_049.{}.fits'.format(x)
orig_file=rho_file

#---- maps file path
b2b_file=lambda x: join(srcdir, 'beta2b//maps/map.unf_{}'.format(x))
a2b_file=lambda x: join(srcdir, 'a2b//maps//map.unf_{}'.format(x))
beta_file=lambda x, y:join(srcdir, 'a2beta//maps/map.unf_gln1_j{}_r{}.unf'.format(x,y))

#---- power spectra file path
nlist_cln=['b2a/cln.unf','beta2b/cln.unf']
nlist_cl=[ 'b2a/meancl.unf','beta2b/meancl.unf']
clx_file=lambda x: join(srcdir, '{}'.format(x))

#----- figure output dir

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
def main():

    #dsx     = 0.396         # pixel size of SDSS detector.
    dsx = 0.05         # pixel size of SDSS detector.
    R  = 2.9918     #vd is velocity dispersion.
    zl = 0.2     #zl is the redshift of the lens galaxy.
    zs = 1.0
    vd = 220    #Velocity Dispersion.
    nnn = 128      #Image dimension
    bsz = dsx*nnn
    nstd = 59

    xx01 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xx02 = np.linspace(-bsz/2.0,bsz/2.0,nnn)+0.5*dsx
    xi2,xi1 = np.meshgrid(xx01,xx02)
    #----------------------------------------------------------------------
    g_amp = 10.0       # peak brightness value
    g_sig = 0.02    # Gaussian "sigma" (i.e., size)
    g_xcen = 0.06   # x position of center (also try (0.0,0.14)
    g_ycen = 0.11    # y position of center
    g_axrat = 1.0   # minor-to-major axis ratio
    g_pa = 0.0      # major-axis position angle (degrees) c.c.w. from x axis
    gpar = np.asarray([g_amp,g_sig,g_xcen,g_ycen,g_axrat,g_pa])
    #----------------------------------------------------------------------
    #g_source = 0.0*xi1
    #g_source = gauss_2d(xi1,xi2,gpar) # modeling source as 2d Gaussian with input parameters.
    #----------------------------------------------------------------------
    xc1 = 0.0       #x coordinate of the center of lens (in units of Einstein radius).
    xc2 = 0.0       #y coordinate of the center of lens (in units of Einstein radius).
    q   = 0.7       #Ellipticity of lens.
    rc  = 0.0       #Core size of lens (in units of Einstein radius).
    re  = re_sv(vd,zl,zs)       #Einstein radius of lens.
    pha = 45.0      #Orintation of lens.
    lpar = np.asarray([xc1,xc2,q,rc,re,pha])
    #----------------------------------------------------------------------
    ai1,ai2,mua = lens_equation_sie(xi1,xi2,lpar)
    yi1 = xi1-ai1
    yi2 = xi2-ai2
    #----------------------------------------------------------------------
    gpar = np.asarray([g_amp,g_sig,g_xcen,g_ycen,g_axrat,g_pa])
    g_limage = gauss_2d(yi1,yi2,gpar)

    #print np.sum(g_limage)
    #pl.figure()
    #pl.contourf(g_limage)
    #pl.colorbar()

    #g_limage = rebin_array_2d.rebin(g_limage,[128,128])
    #-------------------------------------------------------------
    dA = Planck13.comoving_distance(zl).value*1000./(1+zl)
    Re = dA*np.sin(R*np.pi/180./3600.)
    counts  =Brightness(R,vd)
    vpar = np.asarray([counts,Re,xc1,xc2,q,pha])
    #g_lens = deVaucouleurs(xi1,xi2,xc1,xc2,counts,R,1.0-q,pha)
    g_lens = de_vaucouleurs_2d(xi1,xi2,vpar)

    g_lens = g_lens/1e4

    file_psf = "../PSF_and_noise/sdsspsf.fits"
    g_psf = pyfits.getdata(file_psf)-1000.0
    g_psf = g_psf/np.sum(g_psf)
    g_limage = ss.fftconvolve(g_limage+g_lens,g_psf,mode="same")

    g_noise = noise_map(nnn,nnn,nstd,"Gaussian")


    g_limage = g_limage+(g_noise+skycount)/200.0

    pl.figure()
    pl.contourf(g_limage)
    pl.colorbar()

    output_filename = "../output_fits/test.fits"
    pyfits.writeto(output_filename,g_limage,clobber=True)

    pl.show()

    return 0

# 需要導入模塊: from astropy.cosmology import Planck13 [as 別名]
# 或者: from astropy.cosmology.Planck13 import comoving_distance [as 別名]
    twinklesid = agn_twinkles_id[:]

    xs = hdulist[1].data['XSRC'][twinklesid]   # needed from OM10
    ys = hdulist[1].data['YSRC'][twinklesid]
    zl = hdulist[1].data['ZLENS'][twinklesid]

#----------------------------------------------------------------------------
    host_id, host_ra, host_dec, host_mag_norm, host_redshift, host_major_axis, host_minor_axis = np.loadtxt("./twinkles_DESC_SLAC/sprinkled_lens_galaxies_230.txt", comments='#', delimiter=',', converters=None, skiprows=1, usecols=None, unpack=True, ndmin=0)
    host_id = host_id.astype('int')



    # Da_s = p13.angular_diameter_distance(agn_redshift).value
    # Da_l = p13.angular_diameter_distance(zl).value

    Dc_s = p13.comoving_distance(agn_redshift).value
    Dc_l = p13.comoving_distance(zl).value

    idx = agn_image_num == 0

    host_nid = np.array(agn_twinkles_id)[idx]
    lang_nid = np.array(agn_image_num)[idx]
    host_nra = np.array(xs)[idx]
    host_ndec= np.array(ys)[idx]
    host_nmag= np.array(np.random.choice(host_mag_norm, len(agn_id), replace=True)
                        - 5.0*np.log10(Dc_l/Dc_s))[idx]
    host_nred= np.array(agn_redshift)[idx]
    host_nma = np.array(np.random.choice(host_major_axis, len(agn_id)) * Dc_l/Dc_s)[idx]
    host_nmb = np.array(np.random.choice(host_minor_axis, len(agn_id)) * Dc_l/Dc_s)[idx]

    np.savetxt("./twinkles_DESC_SLAC/sprinkled_host_galaxies_230.txt",