LenstronomyMCMC/lenstronomy_test.py at master · Cywtim/LenstronomyMCMC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
# -*- coding: utf-8 -*-
"""
Created on Thu Apr 28 16:45:39 2022

@author: CHENG
"""

from lens_model_class import lens_redshift_difference
import numpy as np
import matplotlib.pyplot as plt
from lenstronomy.LensModel.lens_model import LensModel
from lenstronomy.LensModel.Solver.lens_equation_solver import LensEquationSolver
import time
import copy

#choice of lens mdoel
lens_nie = 'NIE'

lens_model_list = [lens_nie]

#set up model
lensModel = LensModel(lens_model_list=lens_model_list)


#def parameter values of lens models
kwargs_nie1 = {'theta_E':40, 'e1':0.11, 'e2':-0.09,\
              's_scale':1, 'center_x':0, 'center_y':0}
kwargs_nie2 = {'sigma_v':800, 'e1':0.1, 'e2':-0.2,\
              's_scale':0.1, 'center_x':0.10, 'center_y':0.5}


kwargs_lens_list = np.array([kwargs_nie2])


beta = [0.0, 0.0]

cosmo_name = 'w0waCDM'
zl = np.array([0.68])
zs = np.array([1.74])
cosmo_pi = [72,0.3,0.7,-1,0]
s_time = time.time()
lens_model = lens_redshift_difference(lens_model_list)

"""
in_time = time.time()

timedelay_list = lens_model.NIE_timedelay(\
                beta, kwargs_lens_list, cosmo_name, zl, zs, cosmo_pi,\
                  search_window=100,min_distance=3,solver='lenstronomy')
e_time = time.time()

print(timedelay_list.max() - timedelay_list.min())
"""


#print(redshift_difference,',',e_time - s_time)
#redshift_difference = redshift_difference.max() - redshift_difference.min()

wa_list = np.linspace(-1, 1)
redshift_difference_list = []
for i in range(len(wa_list)):
    cosmo_pi[-1] = wa_list[i]
    kwargs_lens = copy.deepcopy(kwargs_lens_list)
    redshift_difference = lens_model.NIE_redshifts(\
                    beta, kwargs_lens, cosmo_name, 0.5, 4, cosmo_pi,\
                      search_window=100,min_distance=3,solver='lenstronomy')
    redshift_difference = redshift_difference.max() - redshift_difference.min()
    redshift_difference_list.append(redshift_difference)


cosmo_pi[-1] = 0
kwargs_lens = copy.deepcopy(kwargs_lens_list)
redshift_difference = lens_model.NIE_redshifts(\
                beta, kwargs_lens, cosmo_name, 0.5, 4, cosmo_pi,\
                  search_window=100,min_distance=3,solver='lenstronomy')
standard = redshift_difference.max() - redshift_difference.min()

plt.figure()
plt.plot(wa_list, np.abs(redshift_difference_list-standard))
plt.xlabel("$w_a$")
plt.ylabel("Redshift difference")


"""
rsd = []
arg_range = np.array([60,130])
for arg in arg_range:
    kwargs_nie1['theta_E'] = arg
    redshift_difference_list = lens_model.lens_redshifts(\
                      beta, kwargs_lens_list, cosmo_name, 1, 4, cosmo_pi,\
                          search_window=300,solver='lenstronomy')

    rsd.append(redshift_difference_list.max() - redshift_difference_list.min())
"""