diff --git a/ChiBenchmark/chiBenchmark/ChiBenchmark b/ChiBenchmark/chiBenchmark/ChiBenchmark
new file mode 100755
index 0000000..614bf34
Binary files /dev/null and b/ChiBenchmark/chiBenchmark/ChiBenchmark differ
diff --git a/ChiBenchmark/chiBenchmark/Makefile b/ChiBenchmark/chiBenchmark/Makefile
new file mode 100644
index 0000000..f9ef30b
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/Makefile
@@ -0,0 +1,50 @@
+PROGRAM_NAME := ChiBenchmark
+#CXX=g++ -lm -ffast-math -ftree-loop-vectorize
+CXX=icpc
+
+
+program_CXX_SRCS := $(wildcard *.cpp)
+program_CXX_OBJS := ${program_CXX_SRCS:.cpp=.o}
+
+#program_CXX_SRCS += $(wildcard ../../*.c) #Find C source files from additonal directories
+program_C_OBJS := ${program_C_SRCS:.c=.o}
+#
+program_CU_SRCS := $(wildcard *.cu)
+#program_CU_SRCS += $(wildcard ../../*.cu) #Find CUDA source files from additional directories
+#program_CU_HEADERS := $(wildcard *.cuh) #Optional: Include .cuh files as dependencies
+#program_CU_HEADERS += $(wildcard ../../*.cuh) #Find .cuh files from additional directories
+program_CU_OBJS := ${program_CU_SRCS:.cu=.cuo}
+#
+program_INCLUDE_DIRS := . /usr/local/cuda/include/ #C++ Include directories
+program_INCLUDE_DIRS += /usr/include/cfitsio/
+program_INCLUDE_DIRS += /users/fgilles/Projects/Libs/cfitsio/include
+program_INCLUDE_DIRS += /users/fgilles/bin/iaca-lin32/include
+program_CU_INCLUDE_DIRS := /home/users/amclaugh/CUB/cub-1.3.2/ #CUDA Include directories
+#
+program_INCLUDE_LIBS := /usr/lib64/ #Include libraries
+program_INCLUDE_LIBS += /users/fgilles/Projects/Libs/cfitsio/lib #Include libraries
+#
+# Compiler flags
+CPPFLAGS += $(foreach includedir,$(program_INCLUDE_DIRS),-I$(includedir))
+CPPFLAGS += $(foreach includelib,$(program_INCLUDE_LIBS),-L$(includelib) -lcfitsio)
+CXXFLAGS += -march=core-avx2 -g -O3 -std=c++0x -Wall -pedantic
+
+
+OBJECTS = $(program_CXX_OBJS) $(program_C_OBJS)
+
+.PHONY: all clean distclean
+#
+all: $(PROGRAM_NAME)
+#
+debug: CXXFLAGS = -g -O0 -std=c++0x -Wall -pedantic -DDEBUG $(EXTRA_FLAGS)
+debug: $(PROGRAM_NAME)
+
+$(PROGRAM_NAME): $(OBJECTS)
+ $(CXX) $(CPPFLAGS) -o $@ $(program_CXX_OBJS) $(program_C_OBJS) $(CUO_O_OBJECTS)
+
+clean:
+ @- $(RM) $(PROGRAM_NAME) $(OBJECTS) *~ *.o *.optrpt
+
+distclean: clean
+
+
diff --git a/ChiBenchmark/chiBenchmark/TestResultParameter.par b/ChiBenchmark/chiBenchmark/TestResultParameter.par
new file mode 100644
index 0000000..dfa5739
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/TestResultParameter.par
@@ -0,0 +1,796 @@
+#The user defines which mode of lenstool he wants to run here
+runmode
+ source 0 sources.cat
+ image 1 theoretical_images_time.txt
+ #inverse 0
+ grid 1 22 1.5
+ poten 1 1000 1.5
+ mass 1 1000 0.4 1.5
+ dpl 1 1000 1.5
+ amplif 1 1000 1.5
+ nbgridcells 100
+ Debug # true activates debug mode
+ end
+frame
+ #dmax 45 #either dmax or x and y have to be declared when a grid is used
+ xmin -50.000
+ xmax 50.000
+ ymin -50.000
+ ymax 50.000
+ end
+cosmology
+ model 1
+ H0 70.000
+ omegaM 0.300
+ omegaX 0.700
+ omegaK 0.000
+ wA 0.000
+ wX -1.000
+ end
+potentiel 1
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 800. #Dispersion Velocity [km/s]
+ rcut 500 #Cut radius (PIEMD distribution)
+ rcore 5 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 2
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 25.000 #X Position [arcsec]
+ y_centre 20.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 140. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 4
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre -16.000 #X Position [arcsec]
+ y_centre 2.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 130. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 5
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre -10.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 110. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 6
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 3.000 #X Position [arcsec]
+ y_centre -10.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 80. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 7
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 17.000 #X Position [arcsec]
+ y_centre -7.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 120. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 8
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 9
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 10
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 11
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 12
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 13
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 14
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 15
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 16
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 17
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 18
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 19
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 20
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 21
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 22
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 23
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 24
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 25
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 26
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 27
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 28
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 29
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 30
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 8
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 9
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 10
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 11
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 12
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 13
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 14
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 15
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 16
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 17
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 18
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 19
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 20
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 21
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 22
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 23
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 24
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 25
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 26
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 27
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 28
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 29
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 30
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 8
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 9
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 10
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 11
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 12
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 13
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 14
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 15
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 16
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 17
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 18
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 19
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 20
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 21
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 22
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 23
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 24
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 25
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 26
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 27
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 28
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 29
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 30
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+limit 1
+ x_centre 1 -1.0 1.0 0.01000
+ y_centre 1 -1.0 1.0 0.010000
+ ellipticity 1 0.05 0.75 0.01
+ angle_pos 0 0. 180.0 0.1
+ rcore 0 0.1 4. 0.100000
+ rcut 0 1.0 20.0 0.10000
+ v_disp 0 200.0 1000.0 0.10000
+ end
+cline
+ nplan 1 1.5
+ dmax 50
+ limitLow 0.2
+ limitHigh 10
+ nbgridcells 1000
+ end
+finish
diff --git a/ChiBenchmark/chiBenchmark/gradient.cpp b/ChiBenchmark/chiBenchmark/gradient.cpp
new file mode 100644
index 0000000..d3f8534
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/gradient.cpp
@@ -0,0 +1,546 @@
+#include <iostream>
+#include <string.h>
+//#include <cuda_runtime.h>
+#include <math.h>
+#include <sys/time.h>
+#include <fstream>
+#include <immintrin.h>
+
+#include "simd_math.h"
+#include "gradient.hpp"
+#include "structure.h"
+//#include "iacaMarks.h"
+
+//#define __INV RCP
+//#define __INV RCP_1NR
+#define __INV RCP_2NR
+//#define __SQRT _mm256_sqrt_pd
+//#define __SQRT SQRT
+//#define __SQRT SQRT_1NR
+#define __SQRT SQRT_2NR
+
+
+
+//
+//
+//
+
+/**** usefull functions for PIEMD profile : see old lenstool ****/
+
+/** I*w,v=0.5 Kassiola & Kovner, 1993 PIEMD, paragraph 4.1
+ *
+ * Global variables used :
+ * - none
+ */
+
+/** @brief This module function calculates profile depended information like the impactparameter b0 and the potential ellipticity epot
+ *
+ * @param lens: mass distribution for which to calculate parameters
+ */
+/*
+void module_readParameters_calculatePotentialparameter(Potential *lens)
+{
+ //std::cout << "module_readParameters_calculatePotentialparameter..." << std::endl;
+ switch (lens->type)
+ {
+
+ case(5): //Elliptical Isothermal Sphere//
+ //impact parameter b0
+ lens->b0 = 4* pi_c2 * lens->vdisp * lens->vdisp ;
+ //ellipticity_potential
+ lens->ellipticity_potential = lens->ellipticity/3 ;
+ break;
+
+ case(8): // PIEMD //
+ //impact parameter b0
+ lens->b0 = 6.*pi_c2 * lens->vdisp * lens->vdisp;
+ //ellipticity_parameter
+ if ( lens->ellipticity == 0. && lens->ellipticity_potential != 0. ){
+ // emass is (a2-b2)/(a2+b2)
+ lens->ellipticity = 2.*lens->ellipticity_potential / (1. + lens->ellipticity_potential * lens->ellipticity_potential);
+ //printf("1 : %f %f \n",lens->ellipticity,lens->ellipticity_potential);
+ }
+ else if ( lens->ellipticity == 0. && lens->ellipticity_potential == 0. ){
+ lens->ellipticity_potential = 0.00001;
+ //printf("2 : %f %f \n",lens->ellipticity,lens->ellipticity_potential);
+ }
+ else{
+ // epot is (a-b)/(a+b)
+ lens->ellipticity_potential = (1. - sqrt(1 - lens->ellipticity * lens->ellipticity)) / lens->ellipticity;
+ //printf("3 : %f %f \n",lens->ellipticity,lens->ellipticity_potential);
+ }
+ break;
+
+ default:
+ std::cout << "ERROR: LENSPARA profil type of clump "<< lens->name << " unknown : "<< lens->type << std::endl;
+ //printf( "ERROR: LENSPARA profil type of clump %s unknown : %d\n",lens->name, lens->type);
+ break;
+ };
+}
+ */
+void module_readParameters_calculatePotentialparameter_SOA(Potential_SOA *lens, int ii)
+{
+ //std::cout << "module_readParameters_calculatePotentialparameter..." << std::endl;
+ switch (lens->type[ii])
+ {
+
+ case(5): /*Elliptical Isothermal Sphere*/
+ //impact parameter b0
+ lens->b0[ii] = 4* pi_c2 * lens->vdisp[ii] * lens->vdisp[ii] ;
+ //ellipticity_potential
+ lens->ellipticity_potential[ii] = lens->ellipticity[ii]/3 ;
+ break;
+
+ case(8): /* PIEMD */
+ //impact parameter b0
+ lens->b0[ii] = 6.*pi_c2 * lens->vdisp[ii] * lens->vdisp[ii];
+ //ellipticity_parameter
+ if ( lens->ellipticity[ii] == 0. && lens->ellipticity_potential[ii] != 0. ){
+ // emass is (a2-b2)/(a2+b2)
+ lens->ellipticity[ii] = 2.*lens->ellipticity_potential[ii] / (1. + lens->ellipticity_potential[ii] * lens->ellipticity_potential[ii]);
+ //printf("1 : %f %f \n",lens->ellipticity[ii],lens->ellipticity_potential[ii]);
+ }
+ else if ( lens->ellipticity[ii] == 0. && lens->ellipticity_potential[ii] == 0. ){
+ lens->ellipticity_potential[ii] = 0.00001;
+ //printf("2 : %f %f \n",lens->ellipticity[ii],lens->ellipticity_potential[ii]);
+ }
+ else{
+ // epot is (a-b)/(a+b)
+ lens->ellipticity_potential[ii] = (1. - sqrt(1 - lens->ellipticity[ii] * lens->ellipticity[ii])) / lens->ellipticity[ii];
+ //printf("3 : %f %f \n",lens->ellipticity[ii],lens->ellipticity_potential[ii]);
+ }
+ break;
+
+ default:
+ std::cout << "ERROR: LENSPARA profil type of clump "<< lens->name << " unknown : "<< lens->type << std::endl;
+ //printf( "ERROR: LENSPARA profil type of clump %s unknown : %d\n",lens->name, lens->type);
+ break;
+ };
+}
+
+
+//
+/**@brief Return the gradient of the projected lens potential for one clump
+ * !!! You have to multiply by dlsds to obtain the true gradient
+ * for the expressions, see the papers : JP Kneib & P Natarajan, Cluster Lenses, The Astronomy and Astrophysics Review (2011) for 1 and 2
+ * and JP Kneib PhD (1993) for 3
+ *
+ * @param pImage point where the result is computed in the lens plane
+ * @param lens mass distribution
+ */
+
+//
+/// Useful functions
+//
+//static
+complex
+piemd_1derivatives_ci05(double x, double y, double eps, double rc)
+{
+ double sqe, cx1, cxro, cyro, rem2;
+ complex zci, znum, zden, zis, zres;
+ double norm;
+ //
+ //std::cout << "piemd_lderivatives" << std::endl;
+ sqe = sqrt(eps);
+ cx1 = (1. - eps) / (1. + eps);
+ cxro = (1. + eps) * (1. + eps);
+ cyro = (1. - eps) * (1. - eps);
+ //
+ rem2 = x * x / cxro + y * y / cyro;
+ /*zci=cpx(0.,-0.5*(1.-eps*eps)/sqe);
+ znum=cpx(cx1*x,(2.*sqe*sqrt(rc*rc+rem2)-y/cx1));
+ zden=cpx(x,(2.*rc*sqe-y));
+ zis=pcpx(zci,lncpx(dcpx(znum,zden)));
+ zres=pcpxflt(zis,b0);*/
+
+ // --> optimized code
+ zci.re = 0;
+ zci.im = -0.5 * (1. - eps * eps) / sqe;
+ znum.re = cx1 * x;
+ znum.im = 2.*sqe * sqrt(rc * rc + rem2) - y / cx1;
+ zden.re = x;
+ zden.im = 2.*rc * sqe - y;
+ norm = zden.re * zden.re + zden.im * zden.im; // zis = znum/zden
+ zis.re = (znum.re * zden.re + znum.im * zden.im) / norm;
+ zis.im = (znum.im * zden.re - znum.re * zden.im) / norm;
+ norm = zis.re;
+ zis.re = log(sqrt(norm * norm + zis.im * zis.im)); // ln(zis) = ln(|zis|)+i.Arg(zis)
+ zis.im = atan2(zis.im, norm);
+ // norm = zis.re;
+ zres.re = zci.re * zis.re - zci.im * zis.im; // Re( zci*ln(zis) )
+ zres.im = zci.im * zis.re + zis.im * zci.re; // Im( zci*ln(zis) )
+ //
+ //zres.re = zis.re*b0;
+ //zres.im = zis.im*b0;
+ //
+ return(zres);
+}
+//
+//// changes the coordinates of point P into a new basis (rotation of angle theta)
+//// y' y x'
+//// * | /
+//// * | / theta
+//// * | /
+//// *|--------->x
+//
+struct point rotateCoordinateSystem(struct point P, double theta)
+{
+ struct point Q;
+
+ Q.x = P.x*cos(theta) + P.y*sin(theta);
+ Q.y = P.y*cos(theta) - P.x*sin(theta);
+
+ return(Q);
+}
+
+__attribute__((noinline))
+static struct point
+grad_halo(const struct point *pImage, const struct Potential *lens)
+{
+ struct point true_coord, true_coord_rotation, result;
+ double R, angular_deviation;
+ complex zis;
+ //
+ //std::cout << "grad_halo..." << lens->type << std::endl;
+ result.x = result.y = 0.;
+ //
+ /*positionning at the potential center*/
+ // Change the origin of the coordinate system to the center of the clump
+ true_coord.x = pImage->x - lens->position.x;
+ true_coord.y = pImage->y - lens->position.y;
+ //
+ switch (lens->type)
+ {
+ case(5): /*Elliptical Isothermal Sphere*/
+ /*rotation of the coordiante axes to match the potential axes*/
+ true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle);
+ //
+ R=sqrt(true_coord_rotation.x*true_coord_rotation.x*(1 - lens->ellipticity/3.) + true_coord_rotation.y*true_coord_rotation.y*(1 + lens->ellipticity/3.)); //ellippot = ellipmass/3
+ result.x = (1 - lens->ellipticity/3.)*lens->b0*true_coord_rotation.x/(R);
+ result.y = (1 + lens->ellipticity/3.)*lens->b0*true_coord_rotation.y/(R);
+ break;
+ case(8): /* PIEMD */
+ /*rotation of the coordiante axes to match the potential axes*/
+ true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle);
+ /*Doing something....*/
+ zis = piemd_1derivatives_ci05(true_coord_rotation.x, true_coord_rotation.y, lens->ellipticity_potential, lens->rcore);
+ //
+ result.x = lens->b0*zis.re;
+ result.y = lens->b0*zis.im;
+ break;
+ default:
+ std::cout << "ERROR: Grad 1 profil type of clump "<< lens->name << " unknown : "<< lens->type << std::endl;
+ break;
+ };
+ return result;
+}
+//
+//
+//
+struct point module_potentialDerivatives_totalGradient(const runmode_param *runmode, const struct point *pImage, const struct Potential *lens)
+{
+ struct point grad, clumpgrad;
+ grad.x = 0;
+ grad.y = 0;
+ std::cout << "nhalos = " << runmode->nhalos << std::endl;
+ for(int i = 0; i < runmode->nhalos; i++)
+ {
+ clumpgrad = grad_halo(pImage, &lens[i]); //compute gradient for each clump separately
+ //nan check
+ //std::cout << clumpgrad.x << " " << clumpgrad.y << std::endl;
+ if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y)
+ {
+ // add the gradients
+ grad.x += clumpgrad.x;
+ grad.y += clumpgrad.y;
+ }
+ }
+ //
+ return(grad);
+}
+
+struct point grad_halo_piemd_SOA(const struct point *pImage, int iterator, const struct Potential_SOA *lens)
+{
+ struct point clumpgrad;
+ clumpgrad.x = 0;
+ clumpgrad.y = 0;
+
+
+ struct point true_coord, true_coord_rotation; //, result;
+ //double R, angular_deviation;
+ complex zis;
+ //
+ //result.x = result.y = 0.;
+ //
+ true_coord.x = pImage->x - lens->position_x[iterator];
+ true_coord.y = pImage->y - lens->position_y[iterator];
+ //
+ // std::cout << "grad_halo..." << lens->type << std::endl;
+ //
+ /*positionning at the potential center*/
+ // Change the origin of the coordinate system to the center of the clump
+ true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle[iterator]);
+ // std::cout << i << ": " << true_coord.x << " " << true_coord.y << " -> "
+ // << true_coord_rotation.x << " " << true_coord_rotation.y << std::endl;
+ //
+ //double sqe, cx1, cxro, cyro, rem2;
+ //
+ double x = true_coord_rotation.x;
+ double y = true_coord_rotation.y;
+ double eps = lens->ellipticity_potential[iterator];
+ double rc = lens->rcore[iterator];
+ //
+ //std::cout << "piemd_lderivatives" << std::endl;
+ //
+ double sqe = sqrt(eps);
+ //
+ double cx1 = (1. - eps)/(1. + eps);
+ double cxro = (1. + eps)*(1. + eps);
+ double cyro = (1. - eps)*(1. - eps);
+ //
+ double rem2 = x*x/cxro + y*y/cyro;
+ //
+ /*zci=cpx(0.,-0.5*(1.-eps*eps)/sqe);
+ znum=cpx(cx1*x,(2.*sqe*sqrt(rc*rc+rem2)-y/cx1));
+ zden=cpx(x,(2.*rc*sqe-y));
+ zis=pcpx(zci,lncpx(dcpx(znum,zden)));
+ zres=pcpxflt(zis,b0);*/
+ // --> optimized code
+ complex zci, znum, zden, zres;
+ double norm;
+ //
+ zci.re = 0;
+ zci.im = -0.5*(1. - eps*eps)/sqe;
+ //
+ znum.re = cx1*x;
+ znum.im = 2.*sqe*sqrt(rc*rc + rem2) - y/cx1;
+ //
+ zden.re = x;
+ zden.im = 2.*rc*sqe - y;
+ norm = (zden.re*zden.re + zden.im*zden.im); // zis = znum/zden
+ //
+ zis.re = (znum.re*zden.re + znum.im*zden.im)/norm;
+ zis.im = (znum.im*zden.re - znum.re*zden.im)/norm;
+ norm = zis.re;
+ zis.re = log(sqrt(norm*norm + zis.im*zis.im)); // ln(zis) = ln(|zis|)+i.Arg(zis)
+ zis.im = atan2(zis.im, norm);
+ // norm = zis.re;
+ zres.re = zci.re*zis.re - zci.im*zis.im; // Re( zci*ln(zis) )
+ zres.im = zci.im*zis.re + zis.im*zci.re; // Im( zci*ln(zis) )
+ //
+ zis.re = zres.re;
+ zis.im = zres.im;
+ //
+ //zres.re = zis.re*b0;
+ //zres.im = zis.im*b0;
+ //
+ //
+ clumpgrad.x = lens->b0[iterator]*zis.re;
+ clumpgrad.y = lens->b0[iterator]*zis.im;
+ //nan check
+ //
+ return(clumpgrad);
+}
+
+struct point grad_halo_sis_SOA(const struct point *pImage, int iterator, const struct Potential_SOA *lens)
+{
+ struct point true_coord, true_coord_rotation, result;
+ double R, angular_deviation;
+ complex zis;
+
+ result.x = result.y = 0.;
+
+ /*positionning at the potential center*/
+ true_coord.x = pImage->x - lens->position_x[iterator]; // Change the origin of the coordinate system to the center of the clump
+ true_coord.y = pImage->y - lens->position_y[iterator];
+
+ /*Elliptical Isothermal Sphere*/
+ /*rotation of the coordiante axes to match the potential axes*/
+ true_coord_rotation = rotateCoordinateSystem(true_coord, lens->ellipticity_angle[iterator]);
+
+ R=sqrt(true_coord_rotation.x*true_coord_rotation.x*(1-lens->ellipticity[iterator]/3.)+true_coord_rotation.y*true_coord_rotation.y*(1+lens->ellipticity[iterator]/3.)); //ellippot = ellipmass/3
+ result.x=(1-lens->ellipticity[iterator]/3.)*lens->b0[iterator]*true_coord_rotation.x/(R);
+ result.y=(1+lens->ellipticity[iterator]/3.)*lens->b0[iterator]*true_coord_rotation.y/(R);
+
+ return result;
+}
+
+struct point module_potentialDerivatives_totalGradient_SOA(const int *Nlens, const struct point *pImage, const struct Potential_SOA *lens)
+{
+ struct point grad, clumpgrad;
+ grad.x=0;
+ grad.y=0;
+
+ //This here could be done with function pointer to better acomodate future ass distributions functions
+ // However I'm unsure of the time of function pointers -> ask gilles
+ //for the moment lens and Nlens is organised the following way : 1. SIS, 2. PIEMD
+
+ //SIS is the first
+ for(int i=0; i<Nlens[0]; i++){
+ clumpgrad=grad_halo_sis_SOA(pImage,i,&lens[0]); //compute gradient for each clump separately
+
+ if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
+ grad.x+=clumpgrad.x;
+ grad.y+=clumpgrad.y;
+ } // add the gradients
+ }
+
+ //PIEMD is the second
+ for(int i=0; i<Nlens[1]; i++){
+ clumpgrad=grad_halo_piemd_SOA(pImage,i,&lens[1]); //compute gradient for each clump separately
+
+ if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
+ grad.x+=clumpgrad.x;
+ grad.y+=clumpgrad.y;
+ } // add the gradients
+ }
+
+ return(grad);
+}
+
+
+struct point grad_halo_piemd_SOA_AVX(const int Nlens, const struct point *pImage, const struct Potential_SOA *lens)
+{
+ struct point grad, clumpgrad;
+ grad.x = 0;
+ grad.y = 0;
+ //std::cout << "optimized gradient: nhalos = " << runmode->nhalos << std::endl;
+ //
+ // smearing the image coordinates on registers
+ __m256d image_x = _mm256_set1_pd(pImage->x);
+ __m256d image_y = _mm256_set1_pd(pImage->y);
+ //
+//#pragma unroll
+ for(int i = 0; i < Nlens; i = i + 4)
+ {
+ //IACA_START;
+ //
+ //std::cout << "starting loop " << std::endl;
+ __m256d two = _mm256_set1_pd(2.);
+ __m256d one = _mm256_set1_pd(1.);
+ __m256d zero = _mm256_set1_pd(0.);
+ __m256d half = _mm256_set1_pd(0.5);
+ __m256d mhalf = _mm256_set1_pd(-0.5);
+ //
+ // 2 loads
+ __m256d true_coord_x = _mm256_sub_pd(image_x, _mm256_loadu_pd(&lens->position_x[i]));
+ __m256d true_coord_y = _mm256_sub_pd(image_y, _mm256_loadu_pd(&lens->position_y[i]));
+ // 2 loafs
+ __m256d rc = _mm256_loadu_pd(&lens->rcore[i]);
+ __m256d b0 = _mm256_loadu_pd(&lens->b0[i]);
+ // 2 loads
+ __m256d eps = _mm256_loadu_pd(&lens->ellipticity_potential[i]);
+ //
+ __m256d one_minus_eps = _mm256_sub_pd(one, eps);
+ __m256d one_plus_eps = _mm256_add_pd(one, eps);
+ __m256d one_plus_eps_rcp = __INV(one_plus_eps);
+ // 1 load
+ __m256d theta = _mm256_loadu_pd(&lens->ellipticity_angle[i]);
+ /*positionning at the potential center*/
+ __m256d cos_theta = _mm256_cos_pd(theta);
+ __m256d sin_theta = _mm256_sin_pd(theta);
+ // rotation: 6 ops
+ __m256d x = _mm256_add_pd(_mm256_mul_pd(true_coord_x, cos_theta), _mm256_mul_pd(true_coord_y, sin_theta));
+ __m256d y = _mm256_sub_pd(_mm256_mul_pd(true_coord_y, cos_theta), _mm256_mul_pd(true_coord_x, sin_theta));
+ //
+ __m256d sqe = __SQRT(eps);
+ //
+ // (1. - eps)/(1. + eps); 3 ops
+ __m256d cx1 = one_minus_eps*one_plus_eps_rcp;
+ // (1. + eps)*(1. + eps); 3 ops
+ __m256d cxro = one_plus_eps*one_plus_eps;
+ // (1. - eps)*(1. - eps); 3 ops
+ __m256d cyro = one_minus_eps*one_minus_eps;
+ //__m256d rem2 = x*x/(cxro) + y*y/(cyro);
+ // ~5 ops
+ __m256d rem2 = x*x*__INV(cxro) + y*y*__INV(cyro);
+ //
+ __m256d zci_re = zero;
+ __m256d zci_im = mhalf*(one - eps*eps)*__INV(sqe); // ~4 ops
+ // 2.*sqe*sqrt(rc*rc + rem2) - y/cx1
+ // 7 ops
+ __m256d znum_re = cx1*x;
+ __m256d znum_im = two*sqe*__SQRT(rc*rc + rem2) - y*__INV(cx1); // ~4 ops
+ //
+ __m256d zden_re = x;
+ __m256d zden_im = _mm256_mul_pd(_mm256_set1_pd(2.), _mm256_mul_pd(rc, sqe));
+ zden_im = _mm256_sub_pd(zden_im, y);
+ //norm = (zden.re*zden.re + zden.im*zden.im); 3 ops
+ __m256d norm = (zden_re*zden_re + zden_im*zden_im);
+ __m256d zis_re = (znum_re*zden_re + znum_im*zden_im)*__INV(norm); // 3 ops
+ __m256d zis_im = (znum_im*zden_re - znum_re*zden_im)*__INV(norm); // 3 ops
+ norm = zis_re;
+ //
+ zis_re = _mm256_log_pd(__SQRT(norm*norm + zis_im*zis_im)); // 3 ops// ln(zis) = ln(|zis|)+i.Arg(zis)
+ zis_re = _mm256_log_pd(__SQRT(norm*norm + zis_im*zis_im)); // 3 ops// ln(zis) = ln(|zis|)+i.Arg(zis)
+ zis_im = _mm256_atan2_pd(zis_im, norm); //
+ //
+ __m256d zres_re = zci_re*zis_re - zci_im*zis_im; // Re( zci*ln(zis) ) 3 ops
+ __m256d zres_im = zci_im*zis_re + zis_im*zci_re; // Im( zci*ln(zis) ) 3 ops
+ //
+ zis_re = zres_re;
+ zis_im = zres_im;
+ //
+ __m256d b0_x, b0_y;
+ //
+ b0_x = _mm256_mul_pd(b0, zis_re); // 1 ops
+ clumpgrad.x = ((double*) &b0_x)[0];
+ clumpgrad.x += ((double*) &b0_x)[1];
+ clumpgrad.x += ((double*) &b0_x)[2];
+ clumpgrad.x += ((double*) &b0_x)[3];
+ //
+ b0_y = _mm256_mul_pd(b0, zis_im); // 1 ops
+ clumpgrad.y = ((double*) &b0_y)[0];
+ clumpgrad.y += ((double*) &b0_y)[1];
+ clumpgrad.y += ((double*) &b0_y)[2];
+ clumpgrad.y += ((double*) &b0_y)[3];
+ //
+ //nan check
+ if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y)
+ {
+ // add the gradients
+ grad.x += clumpgrad.x;
+ grad.y += clumpgrad.y;
+ }
+ }
+ for(int i = Nlens - Nlens%4; i < Nlens; ++i){
+ //std::cout << i << std::endl;
+ clumpgrad=grad_halo_piemd_SOA(pImage,i,&lens[1]); //compute gradient for each clump separately
+ if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y)
+ { //nan check
+ grad.x+=clumpgrad.x;
+ grad.y+=clumpgrad.y;
+ }
+
+
+
+ }
+ //IACA_END;
+ //
+ return(grad);
+}
+
+struct point module_potentialDerivatives_totalGradient_SOA_AVX(const int *Nlens, const struct point *pImage, const struct Potential_SOA *lens)
+{
+ struct point grad, clumpgrad;
+ grad.x = 0;
+ grad.y = 0;
+
+ //SIS is the first
+ for(int i=0; i<Nlens[0]; i++){
+ clumpgrad=grad_halo_sis_SOA(pImage,i,&lens[0]); //compute gradient for each clump separately
+
+ if(clumpgrad.x == clumpgrad.x or clumpgrad.y == clumpgrad.y){ //nan check
+ grad.x+=clumpgrad.x;
+ grad.y+=clumpgrad.y;
+ } // add the gradients
+ }
+
+ clumpgrad = grad_halo_piemd_SOA_AVX(Nlens[1],pImage,&lens[1]);
+ //std::cout << clumpgrad.x << " " << clumpgrad.y << " " << Nlens[1] <<std::endl;
+ grad.x+=clumpgrad.x;
+ grad.y+=clumpgrad.y;
+
+ return grad;
+
+}
diff --git a/ChiBenchmark/chiBenchmark/gradient.hpp b/ChiBenchmark/chiBenchmark/gradient.hpp
new file mode 100644
index 0000000..5f8f892
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/gradient.hpp
@@ -0,0 +1,29 @@
+#pragma once
+#ifndef __GRAD_HPP__
+#define __GRAD_HPP__
+
+#include "structure.h"
+
+/** for both gradient and second derivatives **/
+struct point rotateCoordinateSystem(struct point P, double theta);
+
+/** gradient **/
+struct point module_potentialDerivatives_totalGradient(const runmode_param *runmode, const struct point *pImage, const struct Potential *lens);
+struct point module_potentialDerivatives_totalGradient_SOA(const int *Nlens, const struct point *pImage, const struct Potential_SOA *lens);
+struct point module_potentialDerivatives_totalGradient_SOA_AVX(const int *Nlens, const struct point *pImage, const struct Potential_SOA *lens);
+
+static struct point grad_halo(const struct point *pImage, const struct Potential *lens);
+
+struct point grad_halo_sis_SOA(const struct point *pImage, int iterator,Potential_SOA *lens);
+struct point grad_halo_piemd_SOA(const struct point *pImage, int iterator,Potential_SOA *lens);
+
+struct point grad_halo_sis_SOA_AVX(const struct point *pImage, int iterator,Potential_SOA *lens);
+struct point grad_halo_piemd_SOA_AVX(const struct point *pImage, int iterator,Potential_SOA *lens);
+
+/** PIEMD **/
+complex piemd_1derivatives_ci05(double x, double y, double eps, double rc);
+
+/** Potential **/
+//void module_readParameters_calculatePotentialparameter(Potential *lens);
+void module_readParameters_calculatePotentialparameter_SOA(Potential_SOA *lens, int i);
+#endif
diff --git a/ChiBenchmark/chiBenchmark/main.cpp b/ChiBenchmark/chiBenchmark/main.cpp
new file mode 100644
index 0000000..db3dcd7
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/main.cpp
@@ -0,0 +1,182 @@
+/**
+ * @file main.cpp
+ * @Author Christoph Schaaefer, EPFL (christophernstrerne.schaefer@epfl.ch)
+ * @date October 2016
+ * @brief Benchmark for gradhalo function
+ */
+
+#include <iostream>
+#include <iomanip>
+#include <string.h>
+#include <math.h>
+#include <sys/time.h>
+#include <fstream>
+#include <sys/stat.h>
+//
+#include <mm_malloc.h>
+//
+#include "structure.h"
+#include "timer.h"
+#include "gradient.hpp"
+#include "module_cosmodistances.h"
+#include "module_readParameters.hpp"
+#include "structure.h"
+
+int module_readCheckInput_readInput(int argc, char *argv[])
+{
+/// check if there is a correct number of arguments, and store the name of the input file in infile
+
+ char* infile;
+ struct stat file_stat;
+
+ // If we do not have 3 arguments, stop
+ if ( argc != 3 )
+ {
+ fprintf(stderr, "\nUnexpected number of arguments\n");
+ fprintf(stderr, "\nUSAGE:\n");
+ fprintf(stderr, "lenstool input_file output_directorypath [-n]\n\n");
+ exit(-1);
+ }
+ else if ( argc == 3 )
+ infile=argv[1];
+ std::ifstream ifile(infile,std::ifstream::in); // Open the file
+
+ // check whether the output directory already exists
+ if (stat(argv[2],&file_stat) < 0){
+ mkdir(argv[2], S_IRUSR | S_IWUSR | S_IXUSR | S_IRGRP | S_IWGRP | S_IXGRP | S_IROTH );
+ }
+ else {
+ printf("Error : Directory %s already exists. Specify a non existing directory.\n",argv[2]);
+ exit(-1);
+ }
+
+ // check whether the input file exists. If it could not be opened (ifile = 0), it does not exist
+ if(ifile){
+ ifile.close();
+ }
+ else{
+ printf("The file %s does not exist, please specify a valid file name\n",infile);
+ exit(-1);
+ }
+
+ return 0;
+}
+
+int main(int argc, char *argv[])
+{
+
+
+ // Setting Up the problem
+ //===========================================================================================================
+
+ // This module function reads the terminal input when calling LENSTOOL and checks that it is correct
+ // Otherwise it exits LENSTOOL
+ module_readCheckInput_readInput(argc, argv);
+
+ // This module function reads the cosmology parameters from the parameter file
+ // Input: struct cosmologicalparameters cosmology, parameter file
+ // Output: Initialized cosmology struct
+ cosmo_param cosmology; // Cosmology struct to store the cosmology data from the file
+ std::string inputFile = argv[1]; // Input file
+ module_readParameters_readCosmology(inputFile, cosmology);
+
+ // This module function reads the runmode paragraph and the number of sources, arclets, etc. in the parameter file.
+ // The runmode_param stores the information of what exactly the user wants to do with lenstool.
+ struct runmode_param runmode;
+ module_readParameters_readRunmode(inputFile, &runmode);
+
+ module_readParameters_debug_cosmology(runmode.debug, cosmology);
+ module_readParameters_debug_runmode(runmode.debug, runmode);
+
+
+ //=== Declaring variables
+ struct grid_param frame;
+ struct galaxy images[runmode.nimagestot];
+ struct galaxy sources[runmode.nsets];
+ struct Potential lenses[runmode.nhalos+runmode.npotfile-1];
+ struct cline_param cline;
+ struct potfile_param potfile;
+ struct Potential potfilepotentials[runmode.npotfile];
+ struct potentialoptimization host_potentialoptimization[runmode.nhalos];
+ int nImagesSet[runmode.nsets]; // Contains the number of images in each set of images
+
+ // This module function reads in the potential form and its parameters (e.g. NFW)
+ // Input: input file
+ // Output: Potentials and its parameters
+
+
+ module_readParameters_Potential(inputFile, lenses, runmode.nhalos);
+ module_readParameters_debug_potential(runmode.debug, lenses, runmode.nhalos);
+
+ // This module function reads in the potfiles parameters
+ // Input: input file
+ // Output: Potentials from potfiles and its parameters
+
+ if (runmode.potfile == 1 ){
+ module_readParameters_readpotfiles_param(inputFile, &potfile);
+ module_readParameters_debug_potfileparam(runmode.debug, &potfile);
+ module_readParameters_readpotfiles(&runmode,&potfile,lenses);
+ module_readParameters_debug_potential(runmode.debug, lenses, runmode.nhalos+runmode.npotfile);
+
+ }
+
+
+
+ // This module function reads in the grid form and its parameters (e.g. NFW)
+ // Input: input file
+ // Output: grid and its parameters
+
+
+ module_readParameters_Grid(inputFile, &frame);
+
+
+
+
+ if (runmode.image == 1 or runmode.inverse == 1 or runmode.time > 0){
+
+ // This module function reads in the strong lensing images
+ module_readParameters_readImages(&runmode, images, nImagesSet);
+ //runmode.nsets = runmode.nimagestot;
+ for(int i = 0; i < runmode.nimagestot; ++i){
+
+ images[i].dls = module_cosmodistances_objectObject(lenses[0].z, images[i].redshift, cosmology);
+ images[i].dos = module_cosmodistances_observerObject(images[i].redshift, cosmology);
+ images[i].dr = module_cosmodistances_lensSourceToObserverSource(lenses[0].z, images[i].redshift, cosmology);
+
+ }
+ module_readParameters_debug_image(runmode.debug, images, nImagesSet,runmode.nsets);
+
+ }
+
+ if (runmode.inverse == 1){
+
+ // This module function reads in the potential optimisation limits
+ module_readParameters_limit(inputFile,host_potentialoptimization,runmode.nhalos);
+ module_readParameters_debug_limit(runmode.debug, host_potentialoptimization[0]);
+ }
+
+
+ if (runmode.source == 1){
+ //Initialisation to default values.(Setting sources to z = 1.5 default value)
+ for(int i = 0; i < runmode.nsets; ++i){
+ sources[i].redshift = 1.5;
+ }
+ // This module function reads in the strong lensing sources
+ module_readParameters_readSources(&runmode, sources);
+ //Calculating cosmoratios
+ for(int i = 0; i < runmode.nsets; ++i){
+
+ sources[i].dls = module_cosmodistances_objectObject(lenses[0].z, sources[i].redshift, cosmology);
+ sources[i].dos = module_cosmodistances_observerObject(sources[i].redshift, cosmology);
+ sources[i].dr = module_cosmodistances_lensSourceToObserverSource(lenses[0].z, sources[i].redshift, cosmology);
+ }
+ module_readParameters_debug_source(runmode.debug, sources, runmode.nsets);
+
+ }
+
+
+
+
+
+
+}
diff --git a/ChiBenchmark/chiBenchmark/module_cosmodistances.cpp b/ChiBenchmark/chiBenchmark/module_cosmodistances.cpp
new file mode 100644
index 0000000..9d41a10
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/module_cosmodistances.cpp
@@ -0,0 +1,339 @@
+/**
+* @file module_cosmodistances.cpp
+* @Author Thomas Jalabert, EPFL (me@example.com)
+* @date July 2015
+* @version 0,1
+* @brief Library for the computation of cosmological ratios
+*
+* compute the cosmological ratio of the distances between the lens and the source and the lens and the observer
+*
+*/
+
+
+
+
+/// include header file
+#include <string>
+#include <stdio.h>
+#include <math.h>
+#include <cstring>
+#include <stdlib.h>
+#include "module_cosmodistances.h"
+
+
+
+
+
+// Declare static functions that will only be used in this module
+// The functions are defined further below
+static double module_cosmodistances_cosmo_root(double z,cosmo_param C);
+static double module_cosmodistances_sk(double x, double k);
+static double module_cosmodistances_chi1(double z,cosmo_param C);
+static double module_cosmodistances_chi2(double z1, double z2,cosmo_param C);
+static double module_cosmodistances_chiz(double z,cosmo_param C);
+static double module_cosmodistances_integral_chiz_ab(double a, double b,cosmo_param C);
+
+
+
+// Function defintions
+//==========================================================================================================
+
+/** @brief Calculates ratio distance_(lens-source)/distance_(source)
+* Calculates ratio distance_(lens-source)/distance_(source)
+* @param nsetofimages number of set of images
+* @param nImagesSet set of images
+* @param z_lens redshift of lens
+* @param source sources
+* @param cosmoratio variable where result is stored
+* @param cosmopar cosmological parameter
+*/
+void module_cosmodistances_lensSourceToSource( const int nsetofimages, int nImagesSet[], double z_lens, galaxy source[], double cosmoratio[], cosmo_param cosmopar){
+
+
+ //int imageCounter = 0; // Count the total number of images up to now
+ for(int i=0; i<nsetofimages; i++){
+ //printf("z_lens %f , imag.redshift %f \n" , z_lens,source[0].redshift);
+ cosmoratio[i]=module_cosmodistances_lensSourceToObserverSource(z_lens,source[i].redshift, cosmopar); /// lens efficiency (angular distance)
+ //imageCounter += (nImagesSet[i] - 1); // We add the number of images to skip to get to the next set
+ }
+}
+
+
+
+
+/** @brief Return the angular distance DA(observator(z=0),object(z)) (no unit)
+* Multiply observerObject(z) by c/H0 to get the true value in Mpc.
+* observerObject(z) * c/H0 = DA(0,z) = 1/(1+z) * Sk( integral( 0,z,c*dz/H(z) ) )
+
+* @param z redshift of object
+* @param cosmopar cosmological parameter
+*/
+double module_cosmodistances_observerObject(double z, cosmo_param cosmopar)
+
+{
+ double g;
+
+ if (cosmopar.omegaX == 0.)
+ {
+ g = module_cosmodistances_cosmo_root(z,cosmopar);
+ // Reformulation of the Mattig relation of OL = OK = 0 (De Sitter)
+ return(2.*((1. - cosmopar.omegaM - g)*(1. - g)) / cosmopar.omegaM / cosmopar.omegaM / (1. + z) / (1. + z));
+ }
+ else
+ {
+ if (cosmopar.curvature != 0.)
+ return(module_cosmodistances_sk(module_cosmodistances_chi1(z,cosmopar)*sqrt(fabs(cosmopar.curvature)), cosmopar.curvature)
+ / (1 + z) / sqrt(fabs(cosmopar.curvature)));
+ else
+ return(module_cosmodistances_chi1(z,cosmopar) / (1 + z));
+ }
+}
+
+
+
+
+
+/** @brief Return the angular distance DA(object(z1),object(z2)) divided by c/H0
+* Return the angular distance DA(object(z1),object(z2)) divided by c/H0
+* objectObject(z1,z2) * c/H0 = DA(z1,z2) = 1/(1+z2) * Sk( integral( z1,z2,c*dz/H(z) ) )
+*
+* @param z1 redshift of object 1
+* @param z2 redshift of object 2
+* @param cosmopar cosmological parameter
+*/
+double module_cosmodistances_objectObject(double z1, double z2, cosmo_param cosmopar)
+{
+ double g1, g2;
+
+ if ( z1 >= z2 )
+ return 0.;
+
+ if (cosmopar.omegaX == 0.)
+ {
+ g1 = module_cosmodistances_cosmo_root(z1,cosmopar);
+ g2 = module_cosmodistances_cosmo_root(z2,cosmopar);
+ // Mattig relation for a De Sitter Universe
+ return(2.*((1. - cosmopar.omegaM - g1*g2)*(g1 - g2))
+ / cosmopar.omegaM / cosmopar.omegaM / (1. + z1) / (1. + z2) / (1. + z2));
+ }
+ else
+ {
+ if ( cosmopar.curvature != 0. )
+ return(module_cosmodistances_sk(module_cosmodistances_chi2(z1, z2,cosmopar)*sqrt(fabs(cosmopar.curvature)), cosmopar.curvature)
+ / (1 + z2) / sqrt(fabs(cosmopar.curvature)));
+ else
+ return(module_cosmodistances_chi2(z1, z2,cosmopar) / (1 + z2));
+ }
+}
+
+
+
+
+/****************************************************************/
+
+ /** @brief Return the lens efficacity E=DA(LS) / DA(OS)
+* If zl > zs, return 0.
+*
+* @param zl redshift lens
+* @param zs redshift source
+* @param cosmopar cosmological parameter
+*/
+
+double module_cosmodistances_lensSourceToObserverSource(double zl, double zs, cosmo_param cosmopar)
+{
+ if ( zl >= zs ){
+ printf("*******************\nWarning, a source is between the lens and the observer\n*******************\n");
+ printf("Zl: %f , ZS: %f \n", zl, zs);
+ return (0.);
+ }
+ return (module_cosmodistances_objectObject(zl, zs, cosmopar) / module_cosmodistances_observerObject(zs, cosmopar));
+}
+
+
+
+ /** @brief Calculate square root
+*
+* @param z redshift
+* @param C cosmological parameter
+*/
+
+// Calculate square root
+static double module_cosmodistances_cosmo_root(double z,cosmo_param C)
+{
+ return(sqrt(1. + C.omegaM*z));
+}
+
+
+
+
+// Calculate S_k for cosmology
+ /** @brief Calculate S_k for cosmology
+*
+* @param x
+* @param k curvature
+*/
+static double module_cosmodistances_sk(double x, double k)
+{
+ if (k > 0)
+ return(sin(x));
+ else if (k < 0)
+ return(sinh(x));
+ else
+ return(x);
+}
+
+
+ /** @brief Return 1 / H(z) multiplied by H0
+* chiz(z) / H0 = 1 / H(z) = 1/H0/(1+z)/sqrt( sum(i, Omega_i*(1+z)^(3*(w_i + 1) ) )
+*
+* @param z redshift
+* @param C cosmological parameter
+*/
+
+static double module_cosmodistances_chiz(double z,cosmo_param C)
+{
+ double x;
+ double yy, yyy, y4;
+ double r0, e1, e2, frac;
+
+ x = 1 + z;
+
+
+ switch (C.model) //TV CPL Model
+ {
+ case(1):
+ x = -x * x * C.curvature + x * x * x * C.omegaM + C.omegaX * pow(x, 3 * (1 + C.wX + C.wa)) * exp(-3 * C.wa * z / x);
+ break;
+ case(2): //TV Cardassian (wx is q, wa is n)
+ yy = pow ( (1.+C.curvature)/C.omegaM ,C.wX);
+ yyy = (yy - 1.) * pow( x,3.*C.wX*(C.wa-1.) );
+ y4 = pow(1.+yyy,1./C.wX);
+ x = -x * x * C.curvature + x * x * x * C.omegaM*y4;
+ break;
+ case(3): //TV Interacting DE Model (wa is delta)
+ yy = C.omegaX*pow( x,3.*(1.+C.wX) );
+ yyy = ( C.omegaM/(C.wa+3.*C.wX) ) * ( C.wa*pow(x,3.*(1.+C.wX)) + 3.*C.wX*pow(x,(3.-C.wa)) );
+ x = -x * x * C.curvature + yy +yyy;
+ break;
+ case(4): //TV Holographic Ricci Scale with CPL
+ r0 = C.omegaM/(1.-C.omegaM);
+ e1 = (3./2.)*( (1.+r0+C.wX+4*C.wa)/(1.+r0+3.*C.wa) );
+ e2 = (-1./2.)*( (1.+r0-3.*C.wX)/(1.+r0+3.*C.wa) ) ;
+ frac = ( 1.+r0+3.*C.wa*(x-1.)/x )/(1.+r0) ;
+ yy = pow(x,e1);
+ yyy = pow(frac,e2);
+ x = (yy*yyy)*(yy*yyy);
+ break;
+ default:
+ printf("ERROR: Unknown cosmological model %d\n", C.model);
+ exit(-1);
+
+ }
+
+ if ( x <= 0 )
+ {
+ printf("ERROR : H^2(z)<=0 produced (z,omegaM,omegaX,wX,wa) = (%.3lf,%.3lf,%.3lf,%.3lf,%.3lf)\n",
+ z, C.omegaM, C.omegaX, C.wX, C.wa );
+ exit(-1);
+ }
+ return( 1. / sqrt(x) );
+}
+
+
+
+ /** @brief Return the proper distance Dp(0,z) divided by c/H0
+* chi1(z) * c/H0 = Dp(0,z) = integral( 0, z, c*dz / H(z) )
+*
+* @param z redshift
+* @param C cosmological parameter
+*/
+static double module_cosmodistances_chi1(double z,cosmo_param C)
+{
+ double rez;
+ rez = module_cosmodistances_integral_chiz_ab(0., z,C);
+ return rez;
+}
+
+
+
+ /** @brief Return the proper distance Dp(z1,z2) divided by c/H0
+* chi2(z) * c/H0 = Dp(z1,z2) = integral( z1, z2, c*dz / H(z) )
+*
+* @param zl redshift lens
+* @param zs redshift source
+* @param C cosmological parameter
+*/
+static double module_cosmodistances_chi2(double z1, double z2,cosmo_param C)
+{
+ double rez;
+ rez = module_cosmodistances_integral_chiz_ab(z1, z2,C);
+ return rez;
+}
+
+
+ /** @brief compute the integral of H0/H(z) by trapezes method
+* compute the integral of H0/H(z) by trapezes method
+* @param a,b, cosmologicalparameters cosmopar
+*
+* @param a
+* @param b
+* @param C cosmological parameter
+*/
+static double module_cosmodistances_integral_chiz_ab(double a, double b,cosmo_param C)
+{
+ int i,nit;
+ double res, epsilon=1.e-5; /// accuracy of the integration
+
+ nit=(b-a)/epsilon;
+ res=epsilon*(module_cosmodistances_chiz(a,C)+module_cosmodistances_chiz(b,C))/2.;
+
+ for(i=1;i<nit;i++)
+ {
+ res=res+epsilon*module_cosmodistances_chiz(a+i*epsilon,C);
+ }
+
+ return res;
+}
+
+
+
+
+
+
+ /** @brief Debug function to print the calculated output of the functions
+ * Debug function to print the calculated output of the functions
+*/
+
+// Debug function to print the calculated output of the functions
+
+int module_cosmodistances_debug(int runmode[], double strongLensingRatios_lensSourceToSource[], double weakLensingRatios_lensSourceToSource[], double weakLensing_observerSource[], int numberCleanLens, double cleanlensRatios_lensSourceToSource[], double cleanlens_observerSource[], std::string DEBUG)
+{
+if (strcasecmp(DEBUG.c_str(), "True") == 0) // If we are in debug mode
+{
+
+ if (runmode[0] == 1 or runmode[1] == 1) // If we have strong lensing
+ {
+ printf("DEBUG: Strong lensing cosmo ratios D_LS/D_OS for first 2 sets: %lf, %lf\n\n", strongLensingRatios_lensSourceToSource[0], strongLensingRatios_lensSourceToSource[1]);
+ };
+ if (runmode[2] == 1) // We have weak lensing
+ {
+ printf("DEBUG: Weak lensing D_LS/D_OS for first 2 arclets: %lf, %lf. Weak lensing D_OS for first 2 arclets: %lf, %lf.\n\n", weakLensingRatios_lensSourceToSource[0], weakLensingRatios_lensSourceToSource[1], weakLensing_observerSource[0], weakLensing_observerSource[1]);
+ };
+ if (numberCleanLens > 0) // We have cleanlens mode
+ {
+ printf("DEBUG: Cleanlens D_LS/D_OS for first 2 sources: %lf, %lf. Cleanlens D_OS for first 2 sources: %lf, %lf.\n\n", cleanlensRatios_lensSourceToSource[0], cleanlensRatios_lensSourceToSource[1], cleanlens_observerSource[0], cleanlens_observerSource[1]);
+ };
+
+
+};
+
+return 0;
+}
+
+
+
+
+
+
+
+
diff --git a/ChiBenchmark/chiBenchmark/module_cosmodistances.h b/ChiBenchmark/chiBenchmark/module_cosmodistances.h
new file mode 100644
index 0000000..34d7632
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/module_cosmodistances.h
@@ -0,0 +1,38 @@
+
+/**
+* @file module_cosmodistances.h
+* @date July 2015
+* @version 0,1
+* @brief Header file for cosmoditace calculation module
+*
+* header file
+*
+*/
+
+
+
+
+// Header guard
+#ifndef MODULE_COSMODISTANCES_H
+#define MODULE_COSMODISTANCES_H
+
+
+
+
+// Include
+#include "structure.h"
+#include <string>
+
+
+
+
+
+// Function definitions
+void module_cosmodistances_lensSourceToSource( const int nsetofimages, int nImagesSet[], double z_lens, galaxy image[], double cosmoratio[], cosmo_param cosmopar);
+double module_cosmodistances_observerObject(double z, cosmo_param cosmopar);
+double module_cosmodistances_objectObject(double z1, double z2, cosmo_param cosmopar);
+double module_cosmodistances_lensSourceToObserverSource(double zl, double zs, cosmo_param cosmopar);
+int module_cosmodistances_debug(int runmode[], double strongLensingRatios_lensSourceToSource[], double weakLensingRatios_lensSourceToSource[], double weakLensing_observerSource[], int numberCleanLens, double cleanlensRatios_lensSourceToSource[], double cleanlens_observerSource[], std::string DEBUG );
+
+
+#endif // end header guard
diff --git a/ChiBenchmark/chiBenchmark/module_readParameters.cpp b/ChiBenchmark/chiBenchmark/module_readParameters.cpp
new file mode 100644
index 0000000..593b09e
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/module_readParameters.cpp
@@ -0,0 +1,1655 @@
+/**
+* @file module_readParameters.cpp
+* @Author Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch)
+* @date July 2015
+* @version 0,1
+* @brief read parameters, images, sources and any other relevant information
+*
+* read parameter file
+*
+*/
+
+
+
+// Include
+//===========================================================================================================
+#include <iostream>
+#include <cstring>
+#include <sstream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <fstream>
+#include <cmath>
+#include "structure.h"
+#include "module_readParameters.hpp"
+
+
+
+
+
+//===========================================================================================================
+// Function definitions
+
+/** @brief This module function reads the cosmology parameters from the parameter file
+*
+* This module function reads the cosmology parameters from the parameter file. The given pointer is
+* updated with the results.
+* Read an infile: The structure of the file is as follows:
+ | .
+ | .
+ | .
+ |cosmology
+ | model x <---- this is a value
+ | H0 x
+ | . .
+ | . .
+ | end
+ | .
+ | .
+ | finish
+*
+* @param Cosmology cosmology
+* @param infile path to file
+*/
+
+void module_readParameters_readCosmology(std::string infile, cosmo_param &Cosmology)
+{
+
+ std::string first, second, third, line1, line2;
+
+
+ std::ifstream IN(infile.c_str(),std::ios::in); // open file
+ if ( IN )
+ {
+ std::getline(IN,line1);
+ std::istringstream read1(line1); // create a stream for the line
+ read1 >> first;
+ while(strncmp(first.c_str(), "finish",6) != 0) // read line by line until finish is reached
+ {
+ if ( strncmp(first.c_str(), "cosmolog", 8) == 0){ // if the line contains information about cosmology
+
+ std::getline(IN,line2); // read the line word by word
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+
+ while(strncmp(second.c_str(), "end",3) != 0) // Go ahead until "end" is reached
+ {
+
+ if ( !strcmp(second.c_str(), "model") ) // set model of universe
+ {
+ Cosmology.model=atoi(third.c_str());
+ }
+ else if ( !strcmp(second.c_str(), "H0") ) // set Hubble constant
+ {
+ Cosmology.H0=atof(third.c_str());
+ Cosmology.h = Cosmology.H0 / 50.;
+ }
+ else if ( !strcmp(second.c_str(), "omegaM") || !strcmp(second.c_str(), "omega") ) // set density of matter
+ {
+ Cosmology.omegaM=atof(third.c_str());
+ }
+ else if ( !strcmp(second.c_str(), "omegaX") || !strcmp(second.c_str(), "lambda") ) // set cosmological constant
+ {
+ Cosmology.omegaX=atof(third.c_str());
+ }
+ else if ( !strcmp(second.c_str(), "wX") || !strcmp(second.c_str(), "q") || !strcmp(second.c_str(), "w0") ) // set "q" for Model 2, "wX" for Model 3, "w0" for Model 4
+ {
+ Cosmology.wX=atof(third.c_str());
+ }
+ else if ( !strcmp(second.c_str(), "wa") || !strcmp(second.c_str(), "n") || !strcmp(second.c_str(), "delta") || !strcmp(second.c_str(), "w1") ) // set "n" for Model 2, "delta" for model 3, "w1" for model 4
+ {
+ Cosmology.wa=atof(third.c_str());
+ }
+ else if ( !strcmp(second.c_str(), "omegaK") ) // set universe curvature
+ {
+ Cosmology.curvature=atof(third.c_str());
+ }
+ // read next line
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+
+ }
+
+ // if a flat Universe
+ if ( Cosmology.curvature == 0. )
+ {
+ Cosmology.omegaX = 1 - Cosmology.omegaM;
+ }
+ else
+ Cosmology.curvature = Cosmology.omegaM + Cosmology.omegaX - 1;
+
+ }
+ // read next line
+ std::getline(IN,line1);
+ std::istringstream read1(line1);
+ read1 >> first;
+
+ }
+
+ IN.close();
+
+ }
+ else
+ {
+ fprintf(stderr, "ERROR: file %s not found\n", infile.c_str()); // Exit if file was not found
+ exit(-1);
+ }
+
+}
+
+/** @brief This module function reads the number of sources, arclets, etc. in the parameter file. We need to know this to allocate the memory
+*
+* Function to read the number of multiple images and clumps
+* Check if there is 1 limit for each clump set
+* The program reads the number of potentials and limits defined, and checks whether there are the same number
+* Then it opens the imfile to read the numbers of images and set of images
+* Reads also the size of the multiple images area, the size of a pixel, and the runmode.
+*
+* @param infile path to file
+* @param runmode runmode parameter
+*/
+
+
+void module_readParameters_readRunmode(std::string infile, struct runmode_param *runmode)
+{
+ std::string first, second, third, fourth, fifth, line1, line2;
+ /// Set default values
+ runmode->nbgridcells = 1000;
+ runmode->source = 0;
+ runmode->image = 0;
+ runmode->imagefile = 'Empty';
+ runmode->mass = 0;
+ runmode->mass_gridcells = 1000;
+ runmode->z_mass = 0.4;
+ runmode->z_mass_s = 0.8;
+ runmode->potential = 0;
+ runmode->pot_gridcells = 1000;
+ runmode->z_pot = 0.8;
+ runmode->potfile = 0;
+ runmode->npotfile = 0;
+ runmode->dpl = 0;
+ runmode->dpl_gridcells = 1000;
+ runmode->z_dpl = 0.8;
+ runmode->inverse = 0;
+ runmode->arclet = 0;
+ runmode->debug = 0;
+ runmode->nimagestot = 0;
+ runmode->nsets = 0;
+ runmode->gridcells = 1000;
+ runmode->cline = 0;
+ runmode->time = 0;
+
+ int j=0;
+ std::string imageFile_name;
+ int imageIndex=0;
+ int numberPotentials=0, numberLimits=0;
+
+/*************************** read nhalos, imfile_name, pixel_size, multipleimagesarea_size and runmode from configuration file *********************/
+
+std::ifstream IN(infile.c_str(), std::ios::in);
+ if ( IN )
+ {
+ std::getline(IN,line1);
+ std::istringstream read1(line1); // create a stream for the line
+ read1 >> first; // Read the first word
+ //std::cout<<first;
+ while( strcmp(first.c_str(),"finish") != 0 ) // Continue until we reach finish
+ {
+
+ if ( strncmp(first.c_str(), "runmode", 7) == 0){ // Read in runmode information
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third >> fourth; // Read in 4 words
+ while (strncmp(second.c_str(), "end", 3)) // Read until we reach end
+ {
+ if ( !strcmp(second.c_str(), "nbgridcells") )
+ {
+ sscanf(line2.c_str(), " %*s %d ", &runmode->nbgridcells);
+ //std::cout<< runmode->nbgridcells;
+
+ }
+
+ if ( !strcmp(second.c_str(), "source") )
+ {
+ char filename[FILENAME_SIZE];
+ //std::cout<< line2 << std::endl;
+ sscanf(line2.c_str(), " %*s %d %s ", &runmode->source, &filename);
+ runmode->sourfile = filename;
+ }
+
+ if ( !strcmp(second.c_str(), "image") )
+ {
+ char filename[FILENAME_SIZE];
+ //std::cout<< line2 << std::endl;
+ sscanf(line2.c_str(), " %*s %d %s ", &runmode->image, &filename);
+ runmode->imagefile = filename;
+ //std::cout<< runmode->imagefile << std::endl;
+
+ }
+ if ( !strcmp(second.c_str(), "inverse") )
+ {
+ sscanf(line2.c_str(), " %*s %d ", &runmode->inverse);
+ }
+ if ( !strcmp(second.c_str(), "mass") )
+ {
+ sscanf(line2.c_str(), " %*s %d %d %lf %lf", &runmode->mass, &runmode->mass_gridcells, &runmode->z_mass, &runmode->z_mass_s);
+ }
+ if ( !strcmp(second.c_str(), "poten") )
+ {
+ sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->potential, &runmode->pot_gridcells, &runmode->z_pot);
+ //printf("grid %d, number %d, redshift %f",runmode->potential, runmode->pot_gridcells, runmode->z_pot);
+ }
+ if ( !strcmp(second.c_str(), "dpl") )
+ {
+ sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->dpl, &runmode->dpl_gridcells, &runmode->z_dpl);
+ }
+ if ( !strcmp(second.c_str(), "grid") )
+ {
+ sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->grid, &runmode->gridcells, &runmode->zgrid);
+ }
+ if ( !strcmp(second.c_str(), "amplif") )
+ {
+ sscanf(line2.c_str(), " %*s %d %d %lf", &runmode->amplif, &runmode->amplif_gridcells, &runmode->z_amplif);
+ }
+ if ( !strcmp(second.c_str(), "arclets") )
+ {
+ runmode->arclet = 1; // Not supported yet
+ }
+ if ( !strcmp(second.c_str(), "time") )
+ {
+ sscanf(line2.c_str(), " %*s %d ", &runmode->time);
+
+ }
+ if( !strncmp(second.c_str(), "Debug",5) ) // Read in if we are in debug mode
+ {
+ runmode->debug = 1;
+ }
+
+ // Read the next line
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+ }
+ }
+
+
+ else if (!strncmp(first.c_str(), "potent", 6)) // each time we find a "potential" string in the configuration file, we add an halo
+ {
+ numberPotentials++;
+ }
+ else if (!strncmp(first.c_str(), "limit", 5)) // same for the limit of the potential
+ {
+ numberLimits++;
+ }
+ else if (!strncmp(first.c_str(), "cline", 5))
+ {
+ runmode->cline = 1;
+ }
+ else if (!strncmp(first.c_str(), "potfile", 7))
+ {
+ runmode->potfile = 1;
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third >> fourth; // Read in 4 words
+ while (strncmp(second.c_str(), "end", 3)) // Read until we reach end
+ {
+ if ( !strcmp(second.c_str(), "filein") )
+ {
+ runmode->potfilename = fourth;
+ //std::cout<< line2 << runmode->potfilename;
+ break;
+
+ }
+ // Read the next line
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+ }
+ }
+
+
+ // read the next line
+ std::getline(IN,line1);
+ std::istringstream read1(line1);
+ read1 >> first;
+
+ }
+
+ IN.close();
+
+ //if(numberLimits!=numberPotentials) printf("Warning : Number of clumps different than number of limits in %s\n", infile.c_str()); // must be limits for each halo
+ runmode->nhalos=numberPotentials;
+
+ }
+ else
+ {
+ fprintf(stderr, "ERROR: file %s not found\n", infile.c_str());
+ exit(-1);
+ }
+
+
+
+/******************************** read nset and nimagestot from imfile_name **********************/
+/* the images file must contain
+ 1 x_center y_center a b theta redshift
+ 1 . . . . . .
+ 1
+ 2
+ 2
+ 2
+ 2
+ 2
+ So here we have 3 images in the first set of images, 5 in the second, and 8 images in total
+*/
+
+ /*TODO ImageFile_name is a crude method of getting the image information. Rethink it.*/
+
+ if (runmode->image == 1 or runmode->inverse == 1 or runmode->time >= 1){
+ imageFile_name = runmode->imagefile;
+ //printf("Image to source mode activated\n");
+
+
+ std::ifstream IM(imageFile_name.c_str(), std::ios::in);
+ //printf(" Booo says hello again \n");
+ if ( IM )
+ {
+ while( std::getline(IM,line2) ) // Read every line
+ {
+ j=atoi(line2.c_str());
+ if(j> runmode->nsets){ // If we move on to a new set images, we pick the new index as reference
+ runmode->nsets=j;
+ }
+ imageIndex++;
+ }
+ }
+ else{
+
+ fprintf(stderr, "ERROR: file %s not found\n", imageFile_name.c_str());
+ exit(-1);
+ }
+
+ runmode->nimagestot=imageIndex;
+ IM.close();
+ }
+
+ if (runmode->source == 1){
+ imageFile_name = runmode->sourfile;
+ //printf("Source to image mode activated\n");
+
+ std::ifstream IM(imageFile_name.c_str(), std::ios::in);
+ //printf(" Booo says hello again \n");
+ if ( IM ){
+ int i = 0;
+ while( std::getline(IM,line1) ){ // Read until we reach the end
+ i++;
+ }
+ runmode->nsets = i ;
+ }
+ else{
+
+ fprintf(stderr, "ERROR: file %s not found\n", imageFile_name.c_str());
+ exit(-1);
+ }
+
+ runmode->nimagestot=imageIndex; // Important
+ IM.close();
+
+ }
+
+ if (runmode->potfile == 1){
+ imageFile_name = runmode->potfilename;
+
+ std::ifstream IM(imageFile_name.c_str(), std::ios::in);
+
+ if ( IM ){
+ int i = 0;
+ while( std::getline(IM,line1) ){ // Read until we reach the end
+ if ( line1[0] == '#' ){
+ continue;
+ }
+ i++;
+ }
+ runmode->npotfile = i ;
+ }
+ else{
+
+ fprintf(stderr, "ERROR: file %s not found\n", imageFile_name.c_str());
+ exit(-1);
+ }
+
+ IM.close();
+
+ }
+
+
+printf(" nsets %d , nhalos %d , nimagestot %d npotfile %d \n", runmode->nsets, runmode->nhalos, runmode->nimagestot,runmode->npotfile);
+/*********** read narclets from arclets_filename ***************/
+/* the arclets file must contain
+ id x_center y_center a b theta redshift
+line : 1 . . . . . . .
+ 2
+ 3
+ .
+ .
+ narclets
+*/
+/*
+ std::ifstream IM_arclets(arclets_filename.c_str(), std::ios::in); // Open file
+ if ( IM_arclets )
+ {
+ while( std::getline(IM_arclets,line2) ) // Read every line
+ {
+ arcletIndex++;
+ }
+ narclets=arcletIndex;
+ }
+ else{
+ printf("Error : file %s not found\n",arclets_filename.c_str());
+ exit(-1);
+ }
+ IM_arclets.close();
+
+
+// Read the number of sources in the cleanlens file
+
+ if (cleanlensMode != 0)
+ {
+ std::ifstream IM_cleanlens(cleanlensFile.c_str(), std::ios::in); // Open file
+ if ( IM_cleanlens )
+ {
+ while( std::getline(IM_cleanlens,line2) ) // Read every line
+ {
+ cleanlensIndex++;
+ }
+ numbercleanlens=cleanlensIndex;
+ }
+ else{
+ printf("Error : file %s not found\n",cleanlensFile.c_str());
+ exit(-1);
+ }
+ IM_cleanlens.close();
+ }
+ else
+ {
+ numbercleanlens = 0;
+ }
+ * */
+}
+
+
+/** @brief read the positions, redshifts and numbers of multiple images from the images file
+*
+* This module function reads in the strong lensing images
+* Output: image coordinates, image shapes (semi-minor, semi-major of ellipse and orientation angle), source redshifts, number of images per set
+*
+** the images file must contain
+ 1 x_center y_center major axis minor axis ellipticity angle redshift
+ 1 . . . . . .
+ 1 . . . . . .
+ 2 . . . . . .
+ 2 . . . . . .
+ 2 . . . . . .
+ 2 . . . . . .
+ 2 . . . . . .
+ . . . . . . .
+ . . . . . . .
+ . . . . . . .
+ nsetofimages . . . . . .
+
+ At each line we store in j the index of the set of images to store the redshifts and the number of images of each set
+ At each line we add an images in the position of images' array
+*
+* @param runmode runmode parameter
+* @param image array where images will be stored
+* @param nImagesSet array where the number of images per Set will be stored
+*/
+
+
+/// read the positions, redshifts and numbers of multiple images from the images file
+void module_readParameters_readImages(const struct runmode_param *runmode, struct galaxy image[], int nImagesSet[])
+{
+
+ std::string second, line1;
+ int imageIndex=0;
+ int j=0;
+
+
+
+
+for(int i=0; i<runmode->nsets; i++){
+ nImagesSet[i]=0;
+}
+
+/*********initialisation of nimages array*********/
+for(int i=0; i<runmode->nimagestot; i++){
+
+ /* Initialise here the variables of the images*/
+ image[i].center.x = image[i].center.y = 0;
+ image[i].shape.a = image[i].shape.b = image[i].shape.theta = 0;
+ image[i].redshift = 0;
+
+ }
+
+//printf("imagefile :%d \n", runmode->imagefile.c_str
+
+// Read in images
+ std::ifstream IM(runmode->imagefile.c_str(),std::ios::in);
+ if ( IM )
+ {
+ int i = 0;
+ while( std::getline(IM,line1) ) // Read until we reach the end
+ {
+ // Read in the parameters, * means we read in a parameter but do not store it
+
+ sscanf(line1.c_str(), "%d %lf %lf %lf %lf %lf %lf %lf", &j, &image[i].center.x, &image[i].center.y, &image[i].shape.a, &image[i].shape.b, &image[i].shape.theta, &image[i].redshift, &image[i].mag );
+
+ //Variables
+ nImagesSet[j-1]++; // Increase the counter of the number of system for system with number j-1 by 1
+ imageIndex++;
+ i++;
+ }
+ }
+ else{
+ std::cout << "Error : file " << runmode->imagefile << " not found" << std::endl;
+ exit(-1);
+ }
+
+ IM.close();
+}
+
+/** @brief read the positions, redshifts and numbers of multiple sources from the sources file
+*
+* This module function reads in the strong lensing sources
+* Output: sources coordinates, sources shapes (semi-minor, semi-major of ellipse and orientation angle), source redshifts
+*
+* the images file must contain
+ x_center y_center major axis minor axis ellipticity angle redshift
+ 1 . . . . . .
+ 2 . . . . . .
+ 3 . . . . . .
+
+*
+* @param runmode runmode parameter
+* @param source array where sources will be stored
+*/
+
+
+/// read the positions, redshifts and numbers of multiple images from the images file
+void module_readParameters_readSources(struct runmode_param *runmode, struct galaxy source[])
+{
+
+ std::string second, line1;
+ int j=0;
+
+ //Calculating nset
+ std::ifstream IM(runmode->sourfile.c_str(),std::ios::in);
+ if ( IM ){
+ int i = 0;
+ while( std::getline(IM,line1) ){ // Read until we reach the end
+ i++;
+ }
+ runmode->nsets = i ;
+ }
+ else{
+ std::cout << "Error : file " << runmode->sourfile << " not found" << std::endl;
+ exit(-1);
+ }
+
+ IM.close();
+
+/*********initialisation of nimages array*********/
+for(int i=0; i<runmode->nsets; i++){
+
+ /* Initialise here the variables of the images*/
+ source[i].center.x = source[i].center.y = 0;
+ source[i].shape.a = source[i].shape.b = source[i].shape.theta = 0;
+ source[i].redshift = 0;
+ source[i].mag = 0;
+
+ }
+
+ std::ifstream IM2(runmode->sourfile.c_str(),std::ios::in);
+// Read in images
+ if ( IM2 )
+ {
+ int i = 0;
+ while( std::getline(IM2,line1) ) // Read until we reach the end
+ {
+ // Read in the parameters, * means we read in a parameter but do not store it
+
+ sscanf(line1.c_str(), "%d %lf %lf %lf %lf %lf %lf %lf", &j, &source[i].center.x, &source[i].center.y, &source[i].shape.a, &source[i].shape.b, &source[i].shape.theta, &source[i].redshift , &source[i].mag);
+
+ i++;
+ }
+ }
+ else{
+ std::cout << "Error : file " << runmode->sourfile << " not found" << std::endl;
+ exit(-1);
+ }
+
+ IM2.close();
+}
+
+/** @brief This module function reads the critical and caustic line information
+ *@param infile path to file
+* @param cline cline parameter variable
+*/
+
+void module_readParameters_CriticCaustic(std::string infile, cline_param *cline){
+
+ std::string first, second, third, line1, line2;
+ //Default value initialiasation
+ cline->nplan = 0;
+ for(int i =0; i < NPZMAX; ++i){
+ cline->cz[i] = 0;
+ cline->dos[i] = 0;
+ cline->dls[i] = 0;
+ cline->dlsds[i] = 0;
+ }
+
+ cline->limitLow = 1;
+ cline->dmax = 1;
+ cline->limitHigh = 10;
+ cline->nbgridcells = 1000;
+
+ std::ifstream IN(infile.c_str(), std::ios::in);
+ if ( IN ){
+ while(std::getline(IN,line1)){
+ std::istringstream read1(line1); // create a stream for the line
+ read1 >> first;
+ if ( strncmp(first.c_str(), "cline", 5) == 0){ // Read in runmode information
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third; // Read in 4 words
+ while (strncmp(second.c_str(), "end", 3)) // Read until we reach end
+ {
+ if ( !strcmp(second.c_str(), "nplan") ){
+ sscanf(third.c_str(), "%d", &cline->nplan);
+ if ( cline->nplan > NPZMAX ){
+ cline->nplan = NPZMAX;
+ }
+ int j = 0;
+ while ( read2 >> third )
+ {
+ sscanf(third.c_str(), "%lf", &cline->cz[j]);
+ //printf(" zf %f \n",cline->cz[j]);
+ j++;
+ }
+ }
+
+ if ( !strcmp(second.c_str(), "dmax") )
+ {
+ sscanf(third.c_str(), "%lf", &cline->dmax);
+ cline->xmax = cline->dmax;
+ cline->xmin = -cline->dmax;
+ cline->ymax = cline->dmax;
+ cline->ymin = -cline->dmax;
+ }
+ if ( !strcmp(second.c_str(), "pas") || !strcmp(second.c_str(), "step") || !strcmp(second.c_str(), "limitLow") )
+ {
+ sscanf(third.c_str(), "%lf", &cline->limitLow);
+ }
+ if ( !strcmp(second.c_str(), "limitHigh") )
+ {
+ sscanf(third.c_str(), "%lf", &cline->limitHigh);
+ }
+ if ( !strcmp(second.c_str(), "nbgridcells") )
+ {
+ sscanf(third.c_str(), "%d", &cline->nbgridcells);
+ //std::cout << third << std::endl;
+ //printf("cline->nbgridcells %d", cline->nbgridcells);
+ }
+ // Read the next line
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+ }
+ }
+ } // closes while loop
+
+} // closes if(IN)
+
+
+
+IN.close();
+
+}
+
+
+/** @brief This module function reads the potfile information
+ *@param infile path to file
+* @param potfile_param potfile parameter variable
+*/
+
+void module_readParameters_readpotfiles_param(std::string infile, potfile_param *potfile){
+
+ std::string first, second, third, line1, line2;
+ //Default value potfile initialiasation
+
+ potfile->potid = 0;
+ potfile->ftype = 1;
+ potfile->type = 0;
+ potfile->zlens = 0;
+ potfile->mag0 = 0;
+ potfile->sigma = 0;
+ potfile->core = 0;
+ potfile->corekpc = 0;
+ potfile->ircut = 0;
+ potfile->cut1 = 0;
+ potfile->cut2 = 0;
+ potfile->cutkpc1 = 0;
+ potfile->cutkpc2 = 0;
+ potfile->islope = 0;
+ potfile->slope1 = 0;
+ potfile->slope2 = 0;
+ potfile->npotfile = 0;
+
+ std::ifstream IN(infile.c_str(), std::ios::in);
+ if ( IN ){
+ while(std::getline(IN,line1)){
+ std::istringstream read1(line1); // create a stream for the line
+ read1 >> first;
+ if ( strncmp(first.c_str(), "potfile", 7) == 0){ // Read in potfile information
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third; // Read in 4 words
+ while (strncmp(second.c_str(), "end", 3)) // Read until we reach end
+ {
+
+ if ( !strcmp(second.c_str(), "filein") )
+ {
+ sscanf(line2.c_str(), " %*s %d %s ", &potfile->ftype, potfile->potfile);
+ }
+ else if ( !strcmp(second.c_str(), "type") )
+ {
+ sscanf(third.c_str(), "%d", &potfile->type);
+ }
+ else if ( !strcmp(second.c_str(), "zlens") || !strcmp(second.c_str(), "z_lens"))
+ {
+ sscanf(third.c_str(), "%lf", &potfile->zlens);
+ }
+
+ else if (!strcmp(second.c_str(), "mag0") || !strcmp(second.c_str(), "r200"))
+ {
+ sscanf(third.c_str(), "%lf", &potfile->mag0);
+ }
+ else if (!strcmp(second.c_str(), "select"))
+ {
+ sscanf(third.c_str(), "%d", &potfile->select);
+ }
+ else if (!strcmp(second.c_str(), "core"))
+ {
+ sscanf(third.c_str(), "%lf", &potfile->core);
+ }
+ else if (!strcmp(second.c_str(), "corekpc"))
+ {
+ sscanf(third.c_str(), "%lf", &potfile->corekpc);
+ }
+ else if (!strcmp(second.c_str(), "cut"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ircut, &potfile->cut1, &potfile->cut2);
+ }
+ else if (!strcmp(second.c_str(), "cutkpc"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ircut, &potfile->cutkpc1, &potfile->cutkpc2);
+ }
+ else if (!strcmp(second.c_str(), "slope") || !strcmp(second.c_str(), "m200slope"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->islope, &potfile->slope1, &potfile->slope2);
+ }
+ else if (!strcmp(second.c_str(), "sigma"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->isigma, &potfile->sigma1, &potfile->sigma2);
+ }
+ else if (!strcmp(second.c_str(), "vdslope") || !strcmp(second.c_str(), "c200slope"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ivdslope, &potfile->vdslope1, &potfile->vdslope2);
+ }
+ else if (!strncmp(second.c_str(), "vdscat", 6))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ivdscat, &potfile->vdscat1, &potfile->vdscat2);
+ }
+ else if (!strncmp(second.c_str(), "rcutscat", 8))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ircutscat, &potfile->rcutscat1, &potfile->rcutscat2);
+ }
+ else if (!strcmp(second.c_str(), "a") || !strcmp(second.c_str(), "m200"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ia, &potfile->a1, &potfile->a2);
+ }
+ else if (!strcmp(second.c_str(), "b") || !strcmp(second.c_str(), "c200"))
+ {
+ sscanf(line2.c_str(), " %*s %d%lf%lf", &potfile->ib, &potfile->b1, &potfile->b2);
+ }
+
+ // Read the next line
+ std::getline(IN,line2);
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+ }
+ }
+ } // closes while loop
+
+} // closes if(IN)
+
+
+
+IN.close();
+
+}
+
+/** @brief This module function loads the potential from the potfile into a lens
+ *@param infile path to file
+* @param cline cline parameter variable
+*/
+
+void module_readParameters_readpotfiles(const runmode_param *runmode, potfile_param *potfile, Potential *lens){
+
+ std::string second, line1;
+ double aa,bb;
+
+
+// Read in potentials
+ std::ifstream IM(potfile->potfile,std::ios::in);
+ if ( IM )
+ {
+ int i = runmode->nhalos;
+ while( std::getline(IM,line1) ) // Read until we reach the end
+ {
+
+ // Skip commented lines with #
+ if ( line1[0] == '#' )
+ continue;
+
+ // Default initialisation of clump
+ lens[i].type = potfile->type;
+ lens[i].z = potfile->zlens;
+ lens[i].ellipticity_potential = lens[i].ellipticity = 0.;
+ lens[i].alpha = 0.;
+ lens[i].rcut = 0.;
+ lens[i].rcore = 0.;
+ lens[i].rscale = 0.;
+ lens[i].mag = 0.;
+ lens[i].lum = 0.;
+ lens[i].vdisp = 0.;
+ lens[i].position.x = lens[i].position.y = 0.;
+ lens[i].ellipticity_angle = 0.;
+ lens[i].weight = 0;
+ lens[i].exponent = 0;
+ lens[i].einasto_kappacritic = 0;
+
+
+ // Read a line of the catalog
+ if ( potfile->ftype == 1 )
+ {
+ sscanf( line1.c_str(), "%s%lf%lf%lf%lf%lf%lf%lf",
+ &lens[i].name, &lens[i].position.x, &lens[i].position.y,
+ &aa, &bb, &lens[i].theta, &lens[i].mag, &lens[i].lum);
+ lens[i].ellipticity = (aa*aa-bb*bb)/(aa*aa+bb*bb);
+ if ( lens[i].ellipticity < 0 )
+ {
+ fprintf( stderr, "ERROR: The potfile clump %s has a negative ellipticity.\n", lens[i].name );
+ exit(-1);
+ }
+ //goto NEXT;
+
+ }
+
+
+ //Variables
+ potfile->npotfile++;
+ i++;
+ }
+ }
+ else{
+ std::cout << "Error : file " << potfile->potfile << " not found" << std::endl;
+ exit(-1);
+ }
+
+ IM.close();
+
+}
+
+
+
+/** @brief read the information about arclets
+* !Not used! Will be reworked
+* This module function reads in the arclet images for weak lensing
+* @param Arclet file
+*/
+
+
+void module_readParameters_arclets(std::string arclets_filename, point arclets_position[], ellipse arclets_shape[], double arclets_redshift[])
+{
+ std::string second, line1;
+ int j=0;
+
+
+/******************** read the arclets file *******************/
+/* the arclets file must contain
+ id x_center y_center a b theta redshift
+line : 1 . . . . . . .
+ 2
+ 3
+ .
+ .
+ narclets
+*/
+
+ std::ifstream IM(arclets_filename.c_str(),std::ios::in);
+ if ( IM )
+ {
+ while( std::getline(IM,line1)) // Read every line
+ { // Read in parameters, * means we read the parameter but don't store it (*s)
+ sscanf(line1.c_str(), "%*s %lf %lf %lf %lf %lf %lf", &arclets_position[j].x, &arclets_position[j].y, &arclets_shape[j].a, &arclets_shape[j].b, &arclets_shape[j].theta, &arclets_redshift[j]);
+ j++;
+ }
+ IM.close();
+ }
+ else
+ {
+ printf("Error : file %s not found\n",arclets_filename.c_str());
+ exit(-1);
+ }
+}
+
+
+/** @brief This module function reads in if a parameter will be optimized by the MCMC or stay fixed.
+ *
+ * This module function reads in if a parameter will be optimized by the MCMC or stay fixed.
+ * If it will be optimized, it specifies its minimum and maximum allowed values. Unless declared otherwise by the user, input values are fixed and won't be optimized.
+*
+* read an infile :
+ | .
+ | .
+ |limit
+ | x_center x x x x <--- these values contains : block min max accuracy
+ | y_center x x x x if block=1 this is a free parameter, otherwise not
+ | . . . . . min and max define the extremal value of this parameter
+ | . . . . . accuracy is a criterium of convergence for the MCMC
+ | end
+ | .
+ | .
+ |limit
+ | x_center x x x x
+ | y_center x x x x
+ | . . . . .
+ | . . . . .
+ | end
+ | .
+ | .
+ |finish
+and fills the variables with these values
+* @param infile path to input file
+* @param host_potentialoptimization array where limits will be stored
+* @param nhalos number of mass distributions
+
+*/
+
+void module_readParameters_limit(std::string infile, struct potentialoptimization host_potentialoptimization[], int nhalos )
+{
+ std::string first, second, line1, line2;
+ int i=0;
+
+ double DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */
+
+/*** initialize the block variables to zero (= not to be optimized) ***/
+for(int index=0; index<nhalos; index++)
+{
+ host_potentialoptimization[index].position.x.block=0;
+ host_potentialoptimization[index].position.y.block=0;
+ host_potentialoptimization[index].weight.block=0;
+ host_potentialoptimization[index].ellipticity_angle.block=0;
+ host_potentialoptimization[index].ellipticity.block=0;
+ host_potentialoptimization[index].rcore.block=0;
+ host_potentialoptimization[index].rcut.block=0;
+ host_potentialoptimization[index].rscale.block=0;
+ host_potentialoptimization[index].exponent.block=0;
+ host_potentialoptimization[index].alpha.block=0;
+ host_potentialoptimization[index].einasto_kappacritic.block=0;
+ host_potentialoptimization[index].z.block=0;
+}
+
+
+
+
+// Read in
+std::ifstream IN(infile.c_str(), std::ios::in);
+ if ( IN )
+ {
+ while(std::getline(IN,line1))
+ {
+ std::istringstream read1(line1);
+ read1 >> first;
+
+
+ if (!strncmp(first.c_str(), "limit", 5)) // Read the limits
+ {
+ while(std::getline(IN,line2))
+ {
+ std::istringstream read2(line2);
+ read2 >> second; // Read in 1 word
+ if (strcmp(second.c_str(), "end") == 0) break; // stop reading at "end" and move to next potential limit section
+
+
+ if (!strcmp(second.c_str(), "x_centre") || // Read in for x center
+ !strcmp(second.c_str(), "x_center") )
+ {
+ sscanf(line2.c_str(), "%*s %d %lf %lf %lf", &host_potentialoptimization[i].position.x.block,
+ &host_potentialoptimization[i].position.x.min, &host_potentialoptimization[i].position.x.max, &host_potentialoptimization[i].position.x.sigma);
+ }
+ else if (!strcmp(second.c_str(), "y_centre") || // Read in for y center
+ !strcmp(second.c_str(), "y_center") )
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].position.y.block,
+ &host_potentialoptimization[i].position.y.min, &host_potentialoptimization[i].position.y.max, &host_potentialoptimization[i].position.y.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "v_disp") ) // Read in for ellipticity
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].vdisp.block,
+ &host_potentialoptimization[i].vdisp.min, &host_potentialoptimization[i].vdisp.max, &host_potentialoptimization[i].vdisp.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "ellip_pot") ) // Read in for ellipticity
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].ellipticity_potential.block,
+ &host_potentialoptimization[i].ellipticity_potential.min, &host_potentialoptimization[i].ellipticity_potential.max, &host_potentialoptimization[i].ellipticity_potential.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "ellipticitymass") || !strcmp(second.c_str(), "ellipticity") || !strcmp(second.c_str(), "ellipticite") ) // Read in for ellipticity
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].ellipticity.block,
+ &host_potentialoptimization[i].ellipticity.min, &host_potentialoptimization[i].ellipticity.max, &host_potentialoptimization[i].ellipticity.sigma);
+ }
+ else if (!strcmp(second.c_str(), "ellipticity_angle") || !strcmp(second.c_str(), "angle_pos")) // Read in for ellipticity angle
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].ellipticity_angle.block,
+ &host_potentialoptimization[i].ellipticity_angle.min, &host_potentialoptimization[i].ellipticity_angle.max, &host_potentialoptimization[i].ellipticity_angle.sigma);
+ host_potentialoptimization[i].ellipticity_angle.min *= DTR;
+ host_potentialoptimization[i].ellipticity_angle.max *= DTR;
+ host_potentialoptimization[i].ellipticity_angle.sigma *= DTR;
+ }
+ else if ( !strcmp(second.c_str(), "rcut") || !strcmp(second.c_str(), "cut_radius")) // Read in for r cut
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rcut.block,
+ &host_potentialoptimization[i].rcut.min, &host_potentialoptimization[i].rcut.max, &host_potentialoptimization[i].rcut.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "rcore") || !strcmp(second.c_str(), "core_radius")) // Read in for core radius
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rcore.block,
+ &host_potentialoptimization[i].rcore.min, &host_potentialoptimization[i].rcore.max, &host_potentialoptimization[i].rcore.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "NFW_rs") || // Read in for NFW scale radius
+ !strcmp(second.c_str(), "rscale") )
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].rscale.block,
+ &host_potentialoptimization[i].rscale.min, &host_potentialoptimization[i].rscale.max, &host_potentialoptimization[i].rscale.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "exponent") ) // Read in for exponent
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].exponent.block,
+ &host_potentialoptimization[i].exponent.min, &host_potentialoptimization[i].exponent.max, &host_potentialoptimization[i].exponent.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "alpha") ) // Read in for alpha
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].alpha.block,
+ &host_potentialoptimization[i].alpha.min, &host_potentialoptimization[i].alpha.max, &host_potentialoptimization[i].alpha.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "einasto_kappacritic") ) // Read in for critical kappa
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].einasto_kappacritic.block,
+ &host_potentialoptimization[i].einasto_kappacritic.min, &host_potentialoptimization[i].einasto_kappacritic.max, &host_potentialoptimization[i].einasto_kappacritic.sigma);
+ }
+ else if (!strcmp(second.c_str(), "virial_mass") || // Read in for virial mass
+ !strcmp(second.c_str(), "masse") ||
+ !strcmp(second.c_str(), "m200") ||
+ !strcmp(second.c_str(), "mass"))
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].weight.block,
+ &host_potentialoptimization[i].weight.min, &host_potentialoptimization[i].weight.max, &host_potentialoptimization[i].weight.sigma);
+ }
+ else if ( !strcmp(second.c_str(), "z_lens") ) // Read in for redshift
+ {
+ sscanf(line2.c_str(), "%*s%d%lf%lf%lf", &host_potentialoptimization[i].z.block,
+ &host_potentialoptimization[i].z.min, &host_potentialoptimization[i].z.max, &host_potentialoptimization[i].z.sigma);
+ }
+ } // end of inner while loop
+ i++; // Move to next potential
+ }
+ }
+ }
+IN.close();
+}
+
+
+/** @brief This module function reads in the potential form and its parameters (e.g. NFW).
+*
+* @param infile path to input file
+* @param lens array where mass distributions will be stored
+* @param nhalos number of mass distributions
+*/
+
+
+void module_readParameters_Potential(std::string infile, Potential lens[], int nhalos)
+{
+
+double DTR=acos(-1.)/180.; /* 1 deg in rad = pi/180 */
+ Potential *ilens;
+ std::string first, second, third, line1, line2;
+ int i=0;
+
+std::ifstream IN(infile.c_str(), std::ios::in);
+ if ( IN )
+ {
+ while(std::getline(IN,line1))
+ {
+ std::istringstream read1(line1); // create a stream for the line
+ read1 >> first;
+ if (!strncmp(first.c_str(), "potent", 6)) // Read in potential
+ {
+/***********************************************/
+
+ ilens = &lens[i];
+
+ ilens->position.x = ilens->position.y = 0.;
+ ilens->ellipticity = 0;
+ ilens->ellipticity_potential = 0.;
+ ilens->ellipticity_angle = 0.;
+ ilens->rcut = 0.;
+ ilens->rcore = 0;
+ ilens->weight = 0;
+ ilens->rscale = 0;
+ ilens->exponent = 0;
+ ilens->alpha = 0.;
+ ilens->einasto_kappacritic = 0;
+ ilens->z = 0;
+
+
+ while(std::getline(IN,line2))
+ {
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+ if (strcmp(second.c_str(), "end") == 0) // Move to next potential and initialize it
+ {
+ if ( ilens->z == 0. ) // Check if redshift from current halo was initialized
+ {
+ fprintf(stderr, "ERROR: No redshift defined for potential %d\n", i);
+ exit(-1);
+ }
+
+ break; // Break while loop and move to next potential
+
+ }
+
+ // Read in values
+
+ if ( !strcmp(second.c_str(), "profil") || // Get profile
+ !strcmp(second.c_str(), "profile") )
+ {
+ if(!strcmp(third.c_str(), "PIEMD") ||
+ !strcmp(third.c_str(), "1") )
+ {
+ ilens->type=1;
+ strcpy(ilens->type_name,"PIEMD");//ilens->type_name="PIEMD";
+ }
+ if(!strcmp(third.c_str(), "NFW") ||
+ !strcmp(third.c_str(), "2") )
+ {
+ ilens->type=2;
+ strcpy(ilens->type_name,"NFW");//ilens->type_name="NFW";
+ }
+ if(!strcmp(third.c_str(), "SIES") ||
+ !strcmp(third.c_str(), "3") )
+ {
+ ilens->type=3;
+ strcpy(ilens->type_name,"SIES");//ilens->type_name="SIES";
+ }
+ if(!strncmp(third.c_str(), "point", 5) ||
+ !strcmp(third.c_str(), "4") )
+ {
+ ilens->type=4;
+ strcpy(ilens->type_name,"point");//ilens->type_name="point";
+ }
+ if(!strcmp(third.c_str(), "SIE") ||
+ !strcmp(third.c_str(), "5") )
+ {
+ ilens->type=5;
+ strcpy(ilens->type_name,"SIE");//ilens->type_name="point";
+ }
+ if(!strcmp(third.c_str(), "8") )
+ {
+ ilens->type=8;
+ strcpy(ilens->type_name,"PIEMD1");//ilens->type_name="point";
+ }
+
+ }
+
+ else if (!strcmp(second.c_str(), "name")) // Get name of lens
+ {
+ sscanf(third.c_str(),"%s",ilens->name);
+ }
+
+ else if (!strcmp(second.c_str(), "x_centre") || // Get x center
+ !strcmp(second.c_str(), "x_center") )
+ {
+ ilens->position.x=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "y_centre") || // Get y center
+ !strcmp(second.c_str(), "y_center") )
+ {
+ ilens->position.y=atof(third.c_str());
+ }
+ else if ( !strcmp(second.c_str(), "ellipticitymass") || !strcmp(second.c_str(), "ellipticity") ) // Get ellipticity
+ {
+ ilens->ellipticity=atof(third.c_str());
+ //ilens->ellipticity=ilens->ellipticity/3.;
+ }
+ else if (!strcmp(second.c_str(), "ellipticity_angle") || !strcmp(second.c_str(), "angle_pos")) // Get ellipticity angle
+ {
+ ilens->ellipticity_angle=atof(third.c_str());
+ ilens->ellipticity_angle *= DTR;
+ }
+ else if ( !strcmp(second.c_str(), "rcore") || !strcmp(second.c_str(), "core_radius")) // Get core radius
+ {
+ ilens->rcore=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "rcut") || !strcmp(second.c_str(), "cut_radius")) // Get cut radius
+ {
+ ilens->rcut=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "NFW_rs") || // Get scale radius of NFW
+ !strcmp(second.c_str(), "rscale"))
+ {
+ ilens->rscale=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "exponent") ) // Get exponent
+ {
+ ilens->exponent=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "alpha") ) // Get alpha
+ {
+ ilens->alpha=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "einasto_kappacritic") || // Get critical kappa
+ !strcmp(second.c_str(), "kappacritic"))
+ {
+ ilens->einasto_kappacritic=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "z_lens")) // Get redshift
+ {
+ ilens->z=atof(third.c_str());
+ }
+ else if (!strcmp(second.c_str(), "v_disp")) // Get Dispersion velocity
+ {
+ ilens->vdisp=atof(third.c_str());
+ }
+ else if ( !strncmp(second.c_str(), "virial_mass", 6) || // Get virial mass
+ !strcmp(second.c_str(), "masse") ||
+ !strcmp(second.c_str(), "m200") ||
+ !strcmp(second.c_str(), "mass") )
+ {
+ ilens->weight=atof(third.c_str());
+ }
+
+
+ } // closes inner while loop
+
+ //Calculate parameters like b0, potential ellipticity and anyother parameter depending on the profile
+ module_readParameters_calculatePotentialparameter(ilens);
+
+ i++; // Set counter to next potential
+
+
+
+
+} // closes if loop
+
+} // closes while loop
+
+if(i==0){
+ printf("Parameter “potential” not found in the file %s",infile.c_str());
+ exit(-1);
+ }
+if(i>1){
+ for(int j=1; j<i; j++)
+ {
+ if(lens[j].z!=lens[j-1].z) printf("Warning : Some halos have different redshifts ! They should be in the same plane (same redshift) !\n");
+ }
+ }
+
+} // closes if(IN)
+
+IN.close();
+
+}
+
+/** @brief This module function calculates profile depended information like the impactparameter b0 and the potential ellipticity epot
+ *
+ * @param lens: mass distribution for which to calculate parameters
+*/
+
+void module_readParameters_calculatePotentialparameter(Potential *lens){
+
+ switch (lens->type)
+ {
+
+ case(5): /*Elliptical Isothermal Sphere*/
+ //impact parameter b0
+ lens->b0 = 4* pi_c2 * lens->vdisp * lens->vdisp ;
+ //ellipticity_potential
+ lens->ellipticity_potential = lens->ellipticity/3 ;
+ break;
+
+ case(8): /* PIEMD */
+ //impact parameter b0
+ lens->b0 = 6.*pi_c2 * lens->vdisp * lens->vdisp;
+ //ellipticity_parameter
+ if ( lens->ellipticity == 0. && lens->ellipticity_potential != 0. ){
+ // emass is (a2-b2)/(a2+b2)
+ lens->ellipticity = 2.*lens->ellipticity_potential / (1. + lens->ellipticity_potential * lens->ellipticity_potential);
+ printf("1 : %f %f \n",lens->ellipticity,lens->ellipticity_potential);
+ }
+ else if ( lens->ellipticity == 0. && lens->ellipticity_potential == 0. ){
+ lens->ellipticity_potential = 0.00001;
+ printf("2 : %f %f \n",lens->ellipticity,lens->ellipticity_potential);
+ }
+ else{
+ // epot is (a-b)/(a+b)
+ lens->ellipticity_potential = (1. - sqrt(1 - lens->ellipticity * lens->ellipticity)) / lens->ellipticity;
+ printf("3 : %f %f \n",lens->ellipticity,lens->ellipticity_potential);
+ }
+ break;
+
+ default:
+ printf( "ERROR: LENSPARA profil type of clump %s unknown : %d\n",lens->name, lens->type);
+ break;
+ };
+
+}
+
+
+
+/** @brief This module function reads the grid information
+ *
+* @param infile path to input file
+* @param grid array where grid information will be stored
+*/
+
+
+void module_readParameters_Grid(std::string infile, grid_param *grid)
+{
+
+ std::string first, second, third, line1, line2;
+ double dmax ;
+
+std::ifstream IN(infile.c_str(), std::ios::in);
+ if ( IN )
+ {
+ while(std::getline(IN,line1))
+ {
+ std::istringstream read1(line1); // create a stream for the line
+ read1 >> first;
+ if (!strncmp(first.c_str(), "grille", 7) || !strncmp(first.c_str(), "frame", 5) ) // Read in potential
+ {
+
+ while(std::getline(IN,line2))
+ {
+ std::istringstream read2(line2);
+ read2 >> second >> third;
+ if (strcmp(second.c_str(), "end") == 0) // Move to next potential and initialize it
+ {
+ break; // Break while loop and move to next potential
+ }
+
+ // Read in values
+
+ if (!strcmp(second.c_str(), "dmax"))
+ {
+ sscanf(third.c_str(), "%lf", &dmax);
+ //printf( "\t%s\t%lf\n", second.c_str(), dmax);
+ grid->xmin = -dmax;
+ grid->xmax = dmax;
+ grid->ymin = -dmax;
+ grid->ymax = dmax;
+ grid->rmax = dmax * sqrt(2.);
+ }
+ else if (!strcmp(second.c_str(), "xmin"))
+ {
+ sscanf(third.c_str(), "%lf", &grid->xmin);
+ //printf("\t%s\t%lf\n", second.c_str(), grid->xmin);
+ }
+ else if (!strcmp(second.c_str(), "xmax"))
+ {
+ sscanf(third.c_str(), "%lf", &grid->xmax);
+ //printf("\t%s\t%lf\n", second.c_str(), grid->xmax);
+ }
+ else if (!strcmp(second.c_str(), "ymin"))
+ {
+ sscanf(third.c_str(), "%lf", &grid->ymin);
+ //printf( "\t%s\t%lf\n", second.c_str(), grid->ymin);
+ }
+ else if (!strcmp(second.c_str(), "ymax"))
+ {
+ sscanf(third.c_str(), "%lf", &grid->ymax);
+ //printf( "\t%s\t%lf\n", second.c_str(), grid->ymax);
+ }
+}
+
+} // closes if loop
+
+} // closes while loop
+
+} // closes if(IN)
+
+
+
+IN.close();
+
+}
+
+
+
+/* @brief Get number of sets for cleanlens mode
+ * !!Not used. Will be reworked!!
+ *
+ * This module function reads in the "clean lens" mode sources: These sources are read in and lensed only once assuming a fixed mass model
+// Then we can see the predicted location of the images and compare with their real positions
+*
+* @param clean lens file, number of clean lens images
+* @return coordinates for each image, shape of each image, redshifts, number of sets, number of images per set
+*/
+
+
+// Get number of sets for cleanlens mode
+void module_readParameters_SingleLensingSourcesNumberSets(std::string infile, int &nsetofimages_cleanlens )
+{
+ std::string line1;
+ int setNumber=0;
+ nsetofimages_cleanlens=0;
+
+
+ // Get number of sets of images
+ std::ifstream IM(infile.c_str(),std::ios::in);
+ if ( IM )
+ {
+ while( std::getline(IM,line1))
+ {
+ // Read in parameters
+ sscanf(line1.c_str(), "%d", &setNumber);
+ if(setNumber > nsetofimages_cleanlens) nsetofimages_cleanlens = setNumber;
+ }
+ }
+ IM.close();
+}
+
+/* @brief Read in the position of sources that are just once lensed to see where their images lensed back into the image plane appear and what their shape looks like ("cleanlens", no MCMC)
+ * !!Not used. Will be reworked!!
+ * Read in the position of sources that are just once lensed to see where their images lensed back into the image plane appear and what their shape looks like ("cleanlens", no MCMC)
+*
+* @param clean lens file, number of clean lens images
+* @return coordinates for each image, shape of each image, redshifts, number of sets, number of images per set
+*/
+
+void module_readParameters_SingleLensingSources(std::string infile, point sources[], ellipse sources_shape[], double redshift[], int nimages_cleanlens[], int nsetofimages_cleanlens )
+{
+ std::string line1;
+ int index=0;
+ int setNumber=0;
+
+
+
+
+/*********initialisation of nimages_cleanlens array*********/
+// Get number of images per set
+for(int i=0; i<nsetofimages_cleanlens; i++){
+ nimages_cleanlens[i]=0;
+ }
+
+
+
+
+// Get number of sets of images
+ std::ifstream IM(infile.c_str(),std::ios::in);
+ if ( IM )
+ {
+ while( std::getline(IM,line1))
+ {
+ // Read in parameters
+ sscanf(line1.c_str(), "%d %lf %lf %lf %lf %lf %lf", &setNumber, &sources[index].x, &sources[index].y, &sources_shape[index].a, &sources_shape[index].b, &sources_shape[index].theta, &redshift[index]);
+ nimages_cleanlens[setNumber-1]++; // We increase the number of images for set with number "setNumber-1" by 1 (we start with the 0th set)
+ index++;
+ }
+ }
+ IM.close();
+
+
+
+
+}
+
+
+/** @brief Prints out cosmology
+*/
+
+void module_readParameters_debug_cosmology(int DEBUG, cosmo_param cosmology)
+{
+
+if (DEBUG == 1) // If we are in debug mode
+{
+ printf("\n\n####################### Summary of Input Parameters #######################\n");
+ printf("Cosmology: Cosmology model = %d, omegaM = %lf, omegaX = %lf, curvature = %lf, wX = %lf, wa = %lf, H0 = %lf, h = %lf\n", cosmology.model, cosmology.omegaM, cosmology.omegaX, cosmology.curvature, cosmology.wX, cosmology.wa, cosmology.H0, cosmology.h);
+
+}
+
+}
+
+/** @brief Prints out runmode
+*/
+
+void module_readParameters_debug_runmode(int DEBUG, runmode_param runmode)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+
+ printf("Runmode: nhalos = %d, nsets = %d, nimagestot = %d, source = %d, image = %d, arclet = %d, cline = %d, inverse= %d, DEBUG = %d\n", runmode.nhalos, runmode.nsets, runmode.nimagestot, runmode.source, runmode.image, runmode.arclet, runmode.cline, runmode.inverse, runmode.debug);
+
+}
+
+}
+/** @brief Prints out images
+*/
+void module_readParameters_debug_image(int DEBUG, galaxy image[], int nImagesSet[], int nsets)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+ int index = 0;
+ for ( int i = 0; i < nsets; ++i){
+ for( int j = 0; j < nImagesSet[i]; ++j){
+ printf("Image [%d]: x = %lf, y = %lf, shape: a = %f, b = %f, theta = %lf, redshift = %lf, nImagesSet = %d,\n",index , image[index].center.x, image[index].center.y, image[index].shape.a, image[index].shape.b, image[index].shape.theta, image[index].redshift, nImagesSet[i]);
+ index +=1;
+ //printf( "%d \n", index );
+ }
+ }
+}
+
+}
+/** @brief Prints out source
+*/
+void module_readParameters_debug_source(int DEBUG, galaxy source[], int nsets)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+ for ( int i = 0; i < nsets; ++i){
+ printf("Source[%d]: x = %lf, y = %lf, shape: a = %f, b = %f, theta = %lf, redshift = %lf. \n\n", i, source[i].center.x, source[i].center.y, source[i].shape.a, source[i].shape.b, source[i].shape.theta, source[i].redshift);
+ }
+}
+
+}
+/** @brief Prints out massdistributions
+*/
+void module_readParameters_debug_potential(int DEBUG, Potential lenses[], int nhalos)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+ for ( int i = 0; i < nhalos; ++i){
+ printf("Potential[%d]: x = %f, y = %f, vdisp = %f, type = %d, type_name = %s, name = %s,\n \t ellipticity = %f, ellipticity_pot = %f, ellipticity angle (radians) = %f, rcore = %f, rcut = %f,\n \t rscale = %f, exponent = %f, alpha = %f, einasto kappa critical = %f, z = %f\n", i,lenses[i].position.x, lenses[i].position.y, lenses[i].vdisp, lenses[i].type, lenses[i].type_name, lenses[i].name, lenses[i].ellipticity, lenses[i].ellipticity_potential, lenses[i].ellipticity_angle, lenses[i].rcore, lenses[i].rcut, lenses[i].rscale, lenses[i].exponent, lenses[i].alpha, lenses[i].einasto_kappacritic, lenses[i].z);
+ }
+}
+
+}
+/** @brief Prints out potfile_param
+*/
+void module_readParameters_debug_potfileparam(int DEBUG, potfile_param *potfile)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+
+ printf("Potfile: potid = %d, filename = %s, ftype = %d, type = %d, zlens = %f, mag0 = %f, sigma = %f,\n \t core = %f, corekpc = %f, ircut = %d, cut1 = %f, cut2 = %f,\n \t cutkpc1 = %f, cutkpc2 = %f, islope = %d, slope1 = %f, slope2 = %f\n", potfile->potid,potfile->potfile,potfile->ftype,potfile->type,potfile->zlens,potfile->mag0,potfile->sigma,potfile->core,potfile->corekpc,potfile->ircut,potfile->cut1,potfile->cut2,potfile->cutkpc1,potfile->cutkpc2,potfile->islope,potfile->slope1,potfile->slope2);
+}
+
+}
+/** @brief Prints out cline
+*/
+void module_readParameters_debug_criticcaustic(int DEBUG, cline_param cline)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+
+ printf("Cline: Number of planes = %d ", cline.nplan);
+ for( int i=0; i < cline.nplan; ++i){
+ printf(" z[%d] = %f, ",i ,cline.cz[i]);
+ }
+ printf(" dmax = %f, Low = %f, High= %f, Nb of gridcells %d \n", cline.dmax , cline.limitLow ,cline.limitHigh, cline.nbgridcells);
+
+}
+
+}
+/** @brief Prints out limits
+*/
+void module_readParameters_debug_limit(int DEBUG, struct potentialoptimization host_potentialoptimization)
+{
+if (DEBUG == 1) // If we are in debug mode
+{
+
+ printf("DEBUG: Center.x B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.position.x.block, host_potentialoptimization.position.x.min, host_potentialoptimization.position.x.max, host_potentialoptimization.position.x.sigma);
+ printf("DEBUG: Center.y B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.position.y.block, host_potentialoptimization.position.y.min, host_potentialoptimization.position.y.max, host_potentialoptimization.position.y.sigma);
+ printf("DEBUG: weight.y B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.weight.block, host_potentialoptimization.weight.min, host_potentialoptimization.weight.max, host_potentialoptimization.weight.sigma);
+ printf("DEBUG: ellipticity_angle B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.ellipticity_angle.block, host_potentialoptimization.ellipticity_angle.min, host_potentialoptimization.ellipticity_angle.max, host_potentialoptimization.ellipticity_angle.sigma);
+ printf("DEBUG: ellipticity B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.ellipticity.block, host_potentialoptimization.ellipticity.min, host_potentialoptimization.ellipticity.max, host_potentialoptimization.ellipticity.sigma);
+ printf("DEBUG: ellipticity_potential B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.ellipticity_potential.block, host_potentialoptimization.ellipticity_potential.min, host_potentialoptimization.ellipticity_potential.max, host_potentialoptimization.ellipticity_potential.sigma);
+ printf("DEBUG: rcore B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.rcore.block, host_potentialoptimization.rcore.min, host_potentialoptimization.rcore.max, host_potentialoptimization.rcore.sigma);
+ printf("DEBUG: rcut B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.rcut.block, host_potentialoptimization.rcut.min, host_potentialoptimization.rcut.max, host_potentialoptimization.rcut.sigma);
+ printf("DEBUG: rscale B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.rscale.block, host_potentialoptimization.rscale.min, host_potentialoptimization.rscale.max, host_potentialoptimization.rscale.sigma);
+ printf("DEBUG: exponent B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.exponent.block, host_potentialoptimization.exponent.min, host_potentialoptimization.exponent.max, host_potentialoptimization.exponent.sigma);
+ printf("DEBUG: vdisp B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.vdisp.block, host_potentialoptimization.vdisp.min, host_potentialoptimization.vdisp.max, host_potentialoptimization.vdisp.sigma);
+ printf("DEBUG: alpha B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.alpha.block, host_potentialoptimization.alpha.min, host_potentialoptimization.alpha.max, host_potentialoptimization.alpha.sigma);
+ printf("DEBUG: z B = %d, min = %f, max = %f, sigma = %f \n", host_potentialoptimization.z.block, host_potentialoptimization.z.min, host_potentialoptimization.z.max, host_potentialoptimization.z.sigma);
+
+}
+
+}
+
+
+
+
+
+
diff --git a/ChiBenchmark/chiBenchmark/module_readParameters.hpp b/ChiBenchmark/chiBenchmark/module_readParameters.hpp
new file mode 100644
index 0000000..5f69d57
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/module_readParameters.hpp
@@ -0,0 +1,60 @@
+/**
+* @file module_readParameters.h
+* @Author Thomas Jalabert, EPFL (me@example.com), Christoph Schaaefer, EPFL (christophernstrerne.schaefer@epfl.ch)
+* @date July 2015
+* @version 0,1
+* @brief Header file for readParameter module
+*
+* read parameter file
+*
+*/
+
+// header guard
+#ifndef MODULE_READPARAMETERS_H
+#define MODULE_READPARAMETERS_H
+
+
+
+
+// Include
+#include <iostream>
+//#include "module_readParameters.hpp"
+#include <cstring>
+#include <sstream>
+#include <stdlib.h>
+#include "structure.h"
+
+
+
+
+// Function declarations
+void module_readParameters_readCosmology(std::string infile, cosmo_param &Cosmology);
+void module_readParameters_readRunmode(std::string infile, struct runmode_param *Runmode_param);
+void module_readParameters_readImages(const struct runmode_param *runmode, galaxy image[], int nImagesSet[]);
+void module_readParameters_readSources(struct runmode_param *runmode, struct galaxy source[]);
+void module_readParameters_Grid(std::string infile, grid_param *grid);
+void module_readParameters_CriticCaustic(std::string infile, cline_param *cline);
+void module_readParameters_arclets(std::string arclets_filename, point arclets_position[], ellipse arclets_shape[], double arclets_redshift[]);
+void module_readParameters_limit(std::string infile, struct potentialoptimization host_potentialoptimization[], int nhalos );
+void module_readParameters_Potential(std::string infile, Potential lens[], int nhalos);
+void module_readParameters_calculatePotentialparameter(Potential *ilens);
+void module_readParameters_SingleLensingSourcesNumberSets(std::string infile, int &nsetofimages_cleanlens );
+void module_readParameters_SingleLensingSources(std::string infile, point sources[], ellipse sources_shape[], double redshift[], int nimages_cleanlens[], int nsetofimages_cleanlens );
+void module_readParameters_readpotfiles_param(std::string infile, potfile_param *potfile);
+void module_readParameters_readpotfiles(const runmode_param *runmode, potfile_param *potfile, Potential *lens);
+
+void module_readParameters_debug_cosmology(int DEBUG, cosmo_param cosmology);
+void module_readParameters_debug_runmode(int DEBUG, runmode_param runmode);
+void module_readParameters_debug_image(int DEBUG, galaxy image[],int nImagesSet[],int nsets);
+void module_readParameters_debug_source(int DEBUG, galaxy source[], int nsets);
+void module_readParameters_debug_potential(int DEBUG, Potential potential[], int nhalos);
+void module_readParameters_debug_potfileparam(int DEBUG, potfile_param *potfile);
+void module_readParameters_debug_criticcaustic(int DEBUG, cline_param cline);
+void module_readParameters_debug_limit(int DEBUG, struct potentialoptimization host_potentialoptimization);
+
+
+
+
+#endif // enf of header guard
+
+
diff --git a/ChiBenchmark/chiBenchmark/simd_math.h b/ChiBenchmark/chiBenchmark/simd_math.h
new file mode 100644
index 0000000..f879c72
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/simd_math.h
@@ -0,0 +1,99 @@
+#ifndef SIMD_MATH
+#define SIMD_MATH_
+//
+#include <immintrin.h>
+
+#ifdef __INTEL_COMPILER
+inline __m256d operator + (__m256d a, __m256d b) {return _mm256_add_pd(a, b);}
+inline __m256d operator - (__m256d a, __m256d b) {return _mm256_sub_pd(a, b);}
+inline __m256d operator * (__m256d a, __m256d b) {return _mm256_mul_pd(a, b);}
+inline __m256d operator / (__m256d a, __m256d b) {return _mm256_div_pd(a, b);}
+#endif
+//
+inline __m256d RCP(const __m256d d)
+{
+ const __m128 b = _mm256_cvtpd_ps(d);
+ const __m128 rcp = _mm_rcp_ps (b);
+ __m256d x0 = _mm256_cvtps_pd(rcp);
+ //
+ return x0;
+}
+//
+//
+//
+inline __m256d RCP_1NR(const __m256d d)
+{
+ const __m128 b = _mm256_cvtpd_ps(d);
+ const __m128 rcp = _mm_rcp_ps (b);
+ __m256d x0 = _mm256_cvtps_pd(rcp);
+ //
+ x0 = x0 + x0 - d*x0*x0;
+ // //
+ return x0;
+}
+//
+//
+//
+inline __m256d RCP_2NR(const __m256d d)
+{
+ const __m128 b = _mm256_cvtpd_ps(d);
+ const __m128 rcp = _mm_rcp_ps (b);
+ __m256d x0 = _mm256_cvtps_pd(rcp);
+ //
+ x0 = x0 + x0 - d*x0*x0;
+ x0 = x0 + x0 - d*x0*x0;
+ //
+ return x0;
+}
+
+
+inline __m256d SQRT(const __m256d d)
+{
+ const __m128 b = _mm256_cvtpd_ps(d);
+ const __m128 rcp = _mm_sqrt_ps (b);
+ __m256d x0 = _mm256_cvtps_pd(rcp);
+
+ return x0;
+}
+//
+//
+//
+inline __m256d SQRT_1NR(const __m256d d)
+{
+ const __m128 b = _mm256_cvtpd_ps(d);
+ const __m128 rcp = _mm_sqrt_ps (b);
+ __m256d x0 = _mm256_cvtps_pd(rcp);
+
+ __m256d half = _mm256_set1_pd(0.5);
+ __m256d three = _mm256_set1_pd(3.);
+
+ __m256d a = RCP_1NR(d);
+
+ x0 = half*x0*(three - x0*x0*a);
+
+ return x0;
+}
+//
+//
+//
+inline __m256d SQRT_2NR(const __m256d d)
+{
+ const __m128 b = _mm256_cvtpd_ps(d);
+ const __m128 rcp = _mm_sqrt_ps (b);
+ __m256d x0 = _mm256_cvtps_pd(rcp);
+
+ __m256d half = _mm256_set1_pd(0.5);
+ __m256d three = _mm256_set1_pd(3.);
+
+ __m256d a = RCP_2NR(d);
+
+ x0 = half*x0*(three - x0*x0*a);
+ x0 = half*x0*(three - x0*x0*a);
+
+ return x0;
+}
+
+
+
+
+#endif
diff --git a/ChiBenchmark/chiBenchmark/structure.h b/ChiBenchmark/chiBenchmark/structure.h
new file mode 100644
index 0000000..20d116a
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/structure.h
@@ -0,0 +1,528 @@
+/**
+* @file structure.h
+* @Author Thomas Jalabert, EPFL (me@example.com) , Christoph Schaefer, EPFL (christophernstrerne.schaefer@epfl.ch)
+* @date July 2015
+* @version 0,1
+* @brief Header file to define the used structures (e.g. defined structs)
+*
+* @param configuration file (parameters.par)
+* @return Depends on choice in the configuration file, e.g. least chi2 model
+*/
+
+
+// Header guard
+#ifndef STRUCTURE_H
+#define STRUCTURE_H
+
+
+#include <iostream>
+
+
+
+/*****************************************************************/
+/* */
+/* Constants: Will be migrated to constants.h when there are too many of them*/
+/* */
+/*****************************************************************/
+
+
+// GPU definitions
+#define threadsPerBlock 512 // threads for each set of images
+#define MAXIMPERSOURCE 20 // maximum number of multiple images for one source
+#define MAXIM 200 // maximum total number of images treated
+
+// Dimensions definitions
+#define NPZMAX 9 /* maximum number of critical lines in g_cline struct*/
+//#define CLINESIZE 500 /* maximum number of critical and caustic points for Cline mode. Adjust depending on RAM*/
+#define NPOTFILESIZE 2000 //maximum number of potential in potfiles
+#define DMIN 1e-4 // distance minimale de convergence dans le plan image (in arcsec)
+#define NITMAX 100
+
+
+// gNFW definitions
+#define gNFW_ARRAY_SIZE 1800 // Set the dimension of the gnfw table gNFW_ARRAY_SIZE, must be 1800 for the current table file
+
+// Filename definition
+#define FILENAME_SIZE 50 // size of a filename in .par file
+
+//constants definition
+
+#define pi_c2 7.209970e-06 /* pi en arcsecond/ c^2 = 648000/vol/vol */
+#define cH2piG 0.11585881 /* cH0/2piG en g/cm^2 (H0=50) */
+#define cH4piG 0.057929405 /* cH0/4piG en g/cm^2 (H0=50) */
+#define cH0_4piG 2.7730112e-4 /* cH0/4piG en 10^12 M_Sol/kpc^2 (H0=50) */
+#define d0 29.068701 /* vol/h0*1000/rta -- c/H0 en h-1.kpc/arcsecond (h0=50)*/
+
+#define MCRIT12 .2343165 /* c^3/4piGh0/RTA/RTA in 1e12 M_sol/arcsec^2 (h0=50) */
+
+/*****************************************************************/
+/* */
+/* Types definition */
+/* */
+/*****************************************************************/
+
+/** @brief Point: Structure of 2 coordinates
+ *
+ * @param x: X coordinate
+ * @param y: Y coordinate
+ */
+struct point
+{
+ double x;
+ double y;
+};
+
+/** @brief Complex: Structure of 2 doubles
+ * @param re: Real Part
+ * @param im: Imaginary Part
+ */
+struct complex
+{
+ double re;
+ double im;
+};
+
+/** @brief Segment: Structure of two points
+ */
+struct segment
+{
+ point a;
+ point b;
+};
+
+/** @brief triplet: Structure of three points defining a triangle
+ */
+struct triplet
+{
+ struct point a;
+ struct point b;
+ struct point c;
+};
+
+/** @brief bitriplet: Defines two linked triangles (one in source plane and one in image plane)
+ * @param i: Triangle on image plane
+ * @param s: Triangle on source plane
+ */
+struct bitriplet
+{
+ struct triplet i;
+ struct triplet s;
+};
+
+/** @brief contains the table needed to compute the potential derivatives of general NFW profile
+ */
+typedef struct
+{
+ double alpha_now, x_now, kappa, dpl;
+} gNFWdata;
+
+/** @brief Matrix: 2by2 doubles
+ */
+struct matrix
+{
+ double a;
+ double b;
+ double c;
+ double d;
+};
+
+/** @brief ellipse: for shape computation
+ * @param a: semimajor axis
+ * @param b: semiminor axis
+ * @param theta: shape ellipticity
+ */
+struct ellipse
+{
+ double a;
+ double b;
+ double theta;
+};
+
+/** @brief Storage type for sources, lenses and arclets
+ * @param center: position of the center of galaxy
+ * @param shape: shape of galaxy
+ * @param mag: magnitude
+ * @param redshift: redshift
+ * @param dls: Distance lens-source
+ * @param dos: Distance observer-source
+ * @param dr: dls/dos
+ */
+
+struct galaxy
+{
+ //char name[IDSIZE];
+ struct point center;
+ struct ellipse shape;
+ double mag;
+ double redshift;
+ double dls;
+ double dos;
+ double dr;
+};
+
+/** @brief Contains the information for optimising a parameter in the inverse mode
+ * @param block: blockorfree variable (whether a parameter is blocked or free for the mcmc algorithm)
+ * @param min: lower optimisation value
+ * @param max: upper optimisation value
+ * @param sigma: optimisation step (MIGHT NOT BE TAKEN INTO ACCOUNT)
+ */
+struct optimize_block
+{
+ int block;
+ double min;
+ double max;
+ double sigma;
+};
+/** @brief two optimize_block to simulate a point
+ */
+struct optimize_point // blockorfree for the point structure
+{
+ struct optimize_block x;
+ struct optimize_block y;
+};
+
+/** @brief Contains the information for optimising the potential in the inverse mode
+ * @param position: position of the center of the halo
+ * @param weight: weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW)
+ * @param b0: Impact parameter
+ * @param ellipticity_angle: orientation of the clump
+ * @param ellipticity: Mass ellipticity
+ * @param ellipticity_potential: Potential ellipticity
+ * @param rcore: PIEMD specific value
+ * @param rcut: PIEMD specific value
+ * @param rscale: scale radius for NFW, Einasto
+ * @param exponent: exponent for Einasto
+ * @param vdisp: Dispersion velocity
+ * @param alpha: exponent for general NFW
+ * @param einasto_kappacritic: critical kappa for Einasto profile
+ * @param z: redshift
+ */
+struct potentialoptimization // block or free variable for the MCMC for the potential
+{
+ struct optimize_point position;
+ struct optimize_block weight;
+ struct optimize_block b0;
+ struct optimize_block ellipticity_angle;
+ struct optimize_block ellipticity;
+ struct optimize_block ellipticity_potential;
+ struct optimize_block rcore;
+ struct optimize_block rcut;
+ struct optimize_block rscale;
+ struct optimize_block exponent;
+ struct optimize_block vdisp;
+ struct optimize_block alpha;
+ struct optimize_block einasto_kappacritic;
+ struct optimize_block z;
+
+
+};
+
+/** @brief Contains the information of the lens potential
+ * @param type: 1=PIEMD , 2=NFW, 3=SIES, 4=point mass, 5=SIS, 8=PIEMD
+ * @param type_name: IEMD, NFW, SIES, point
+ * @param name: name of the clump (e.g. name of the galaxy) : not compulsory
+ * @param position: position of the center of the halo
+ * @param weight: weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW)
+ * @param b0: Impact parameter
+ * @param ellipticity_angle:
+ * @param ellipticity: Mass ellipticity
+ * @param ellipticity_potential: Potential ellipticity
+ * @param rcore: PIEMD specific value
+ * @param rcut: PIEMD specific value
+ * @param rscale: scale radius for NFW, Einasto
+ * @param exponent: exponent for Einasto
+ * @param vdisp: Dispersion velocity
+ * @param alpha: exponent for general NFW
+ * @param einasto_kappacritic: critical kappa for Einasto profile
+ * @param z: redshift
+ */
+
+struct Potential_SOA
+{
+ int* type; // 1=PIEMD ; 2=NFW; 3=SIES, 4=point mass
+ char type_name[10]; // PIEMD, NFW, SIES, point
+ char name[20]; // name of the clump (e.g. name of the galaxy) : not compulsory
+ //struct point position; // position of the center of the halo
+ double* position_x; // position of the center of the halo
+ double* position_y; // position of the center of the halo
+ double* weight; // weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW)
+ double* b0; // Impact parameter
+ double* vdisp; //Dispersion velocity
+ double* ellipticity_angle; // orientation of the clump
+ double* ellipticity; // ellipticity of the mass distribition
+ double* ellipticity_potential; //ellipticity of the potential
+ double* rcore; // core radius
+ double* rcut; // cut radius
+ double* rscale; // scale radius for NFW, Einasto
+ double* exponent; // exponent for Einasto
+ double* alpha; // exponent for general NFW
+ double* einasto_kappacritic; // critical kappa for Einasto profile
+ double* z; // redshift of the clump
+ double* mag; //magnitude
+ double* lum; //luminosity
+ double* theta; //theta
+ double* sigma; // sigma
+};
+
+struct Potential
+{
+ int type; // 1=PIEMD ; 2=NFW; 3=SIES, 4=point mass
+ char type_name[10]; // PIEMD, NFW, SIES, point
+ char name[20]; // name of the clump (e.g. name of the galaxy) : not compulsory
+ struct point position; // position of the center of the halo
+ double weight; // weight of the clump (the projected mass sigma0 for PIEMD, the density rhoscale for NFW)
+ double b0; // Impact parameter
+ double vdisp; //Dispersion velocity
+ double ellipticity_angle; // orientation of the clump
+ double ellipticity; // ellipticity of the mass distribition
+ double ellipticity_potential; //ellipticity of the potential
+ double rcore; // core radius
+ double rcut; // cut radius
+ double rscale; // scale radius for NFW, Einasto
+ double exponent; // exponent for Einasto
+ double alpha; // exponent for general NFW
+ double einasto_kappacritic; // critical kappa for Einasto profile
+ double z; // redshift of the clump
+ double mag; //magnitude
+ double lum; //luminosity
+ double theta; //theta
+ double sigma; // sigma
+};
+
+/*****************************************************************/
+/* */
+/* Control structure definition */
+/* */
+/*****************************************************************/
+
+/** @brief Control structure for runmode parameters
+ *
+ * Default values are set in module_readParameters_readRunmode
+ *
+ * @param nbgridcells: Number of grid cells
+ * @param source: flag for simple source to image conversion
+ * @param sourfile: file name for source information
+ * @param image: flag for simple image to source conversion
+ * @param imafile: file name for image information
+ * @param mass: flag for mass fitsfile
+ * @param massgridcells: Nb of cell for fitsfile
+ * @param z_mass: redshift for which to be computed
+ * @param z_mass_s: redshift of source for which effect of mass will be computed
+ * @param potential: flag for potential fitsfile
+ * @param potgridcells: Nb of cell for fitsfile
+ * @param z_pot: redshift for which to be computed
+ * @param dpl: flag for displacement fitsfile
+ * @param dplgridcells: Nb of cell for fitsfile
+ * @param z_dpl: redshift for which to be computed
+ * @param inverse: flag for inversion mode (MCMC etc,)
+ * @param arclet: flag for arclet mode
+ * @param debug: flag for debug mode
+ * @param nimagestotal: total number of lensed images in file
+ * @param nsets: number of sources attributed to the lensed images
+ * @param nhalo: Number of halos
+ * @param grid: 0 for automatic grid (not implemented), 1 for grid infor applying on source plane, 2 for grid information applying on image plane
+ * @param nbgridcells: Number of grid cells
+ * @param zgrid: redshift of grid
+ * @param cline: flag for critical and caustic line calculation
+ */
+
+struct runmode_param
+{
+ int nbgridcells;
+ //Source Mode
+ int source;
+ std::string sourfile;
+ int nsets;
+ //Image Mode
+ int image;
+ std::string imagefile;
+ int nimagestot;
+ //Mass Mode
+ int mass;
+ int mass_gridcells;
+ double z_mass;
+ double z_mass_s;
+ //Potential Mode
+ int potential;
+ int pot_gridcells;
+ double z_pot;
+ int nhalos;
+ //Potfile Mode
+ int potfile;
+ int npotfile;
+ std::string potfilename;
+ //displacement Mode
+ int dpl;
+ int dpl_gridcells;
+ double z_dpl;
+ //Inverse Mode
+ int inverse;
+ //Arclet Mode
+ int arclet;
+ //Debug Mode
+ int debug;
+ //Grid Mode
+ int grid;
+ int gridcells;
+ double zgrid;
+ //Critic and caustic mode
+ int cline;
+ //Amplification Mode
+ int amplif;
+ int amplif_gridcells;
+ double z_amplif;
+ //Time/Benchmark mode
+ int time;
+};
+
+/** @brief Not used yet
+ *
+ */
+struct image_param
+{
+
+};
+
+/** @brief Not used yet
+ *
+ */
+struct source_param
+{
+
+};
+
+/** @brief Contains Grid information
+ */
+
+struct grid_param
+{
+ double xmin;
+ double xmax;
+ double ymin;
+ double ymax;
+ double lmin;
+ double lmax;
+ double rmax;
+};
+
+/** @brief Control structure for cosmological parameters
+ *
+ * @param model: Cosmological model
+ * @param omegaM:
+ * @param omegaX:
+ * @param curvature: curvature parameter
+ * @param wX:
+ * @param wa:
+ * @param H0: Hubble constant
+ * @param h: H0/50
+ */
+
+struct cosmo_param
+{
+ int model;
+ double omegaM;
+ double omegaX;
+ double curvature;
+ double wX;
+ double wa;
+ double H0;
+ double h;
+};
+
+/** @brief Control structure for potfile parameters
+ *
+ * @param potid: 1: pot P, 2: pot Q
+ @param ftype:
+ @param potfile[FILENAME_SIZE];
+ @param type;
+ @param zlens;
+ @param core;CCC
+ @param corekpc;
+ @param mag0;
+ @param select;
+ @param ircut;
+ @param cut, cut1, cut2;
+ @param cutkpc1, cutkpc2;
+ @param isigma;
+ @param sigma, sigma1, sigma2;
+ @param islope;
+ @param slope, slope1, slope2;
+ @param ivdslope;
+ @param vdslope, vdslope1, vdslope2;
+ @param ivdscat;
+ @param vdscat, vdscat1, vdscat2;
+ @param ircutscat;
+ @param rcutscat, rcutscat1, rcutscat2;
+ @param ia; // scaling relation of msm200
+ @param a, a1, a2;
+ @param ib; // scaling relation of msm200
+ @param b, b1, b2;
+
+ */
+
+struct potfile_param
+{
+ int potid; // 1: pot P, 2: pot Q
+ int ftype;
+ char potfile[FILENAME_SIZE];
+ int type;
+ double zlens;
+ double core;
+ double corekpc;
+ double mag0;
+ int select;
+ int ircut;
+ double cut, cut1, cut2;
+ double cutkpc1, cutkpc2;
+ int isigma;
+ double sigma, sigma1, sigma2;
+ int islope;
+ double slope, slope1, slope2;
+ int ivdslope;
+ double vdslope, vdslope1, vdslope2;
+ int ivdscat;
+ double vdscat, vdscat1, vdscat2;
+ int ircutscat;
+ double rcutscat, rcutscat1, rcutscat2;
+ int ia; // scaling relation of msm200
+ double a, a1, a2;
+ int ib; // scaling relation of msm200
+ double b, b1, b2;
+ int npotfile;
+};
+
+/** @brief Control structure for caustic and critical lines
+ *
+ * @param nplan: number of sourceplanes for which the caustic lines will be computed
+ * @param cz: redshift values array for respective sourceplanes
+ * @param dos: distcosmo1 to redshift z
+ * @param dls: distcosmo2 between lens[0] and z
+ * @param dlsds: ratio of dl0s/dos
+ * @param limitLow: minimum size of the squares in MARCHINGSQUARES
+ * @param dmax: Size of search area
+ * @param xmin:
+ * @param xmax:
+ * @param ymin:
+ * @param ymax:
+ * @param limithigh: maximum size of the squares in MARCHINGSQUARES algorithm
+ * @param nbgridcells: nbgridcells * nbgridcells = number of pixels for critical line computation
+*/
+
+struct cline_param
+{
+ int nplan;
+ double cz[NPZMAX];
+ double dos[NPZMAX]; // distcosmo1 to redshift z
+ double dls[NPZMAX]; // distcosmo2 between lens[0] and z
+ double dlsds[NPZMAX]; // ratio of dl0s/dos
+ double limitLow; // minimum size of the squares in MARCHINGSQUARES or initial step size in SNAKE
+ double dmax;
+ double xmin;
+ double xmax;
+ double ymin;
+ double ymax;
+ double limitHigh; // maximum size of the squares in MARCHINGSQUARES algorithm
+ int nbgridcells; // nbgridcells * nbgridcells = number of pixels for critical line computation
+};
+
+#endif // if STRUCTURE_H
diff --git a/ChiBenchmark/chiBenchmark/theoretical_images_time.txt b/ChiBenchmark/chiBenchmark/theoretical_images_time.txt
new file mode 100644
index 0000000..f6e95f8
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/theoretical_images_time.txt
@@ -0,0 +1,252 @@
+1 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+1 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+1 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+2 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+2 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+2 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+3 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+3 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+3 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+4 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+4 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+4 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+5 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+5 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+5 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+6 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+6 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+6 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+7 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+8 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+8 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+8 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+9 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+9 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+9 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+10 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+10 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+10 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+11 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+11 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+11 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+12 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+12 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+12 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+13 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+13 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+13 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+14 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+14 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+14 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+15 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+15 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+15 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+16 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+16 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+16 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+17 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+18 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+18 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+18 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+19 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+19 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+19 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+20 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+20 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+20 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+21 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+21 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+21 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+22 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+22 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+22 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+23 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+23 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+23 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+24 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+24 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+24 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+25 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+25 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+25 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+26 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+26 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+26 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+27 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+28 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+28 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+28 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+29 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+29 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+29 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+30 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+30 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+30 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+31 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+31 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+31 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+32 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+32 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+32 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+33 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+33 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+33 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+34 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+34 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+34 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+35 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+35 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+35 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+36 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+36 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+36 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+37 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+38 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+38 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+38 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+39 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+39 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+39 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+40 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+40 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+40 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+41 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+41 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+41 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+42 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+42 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+42 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+43 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+43 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+43 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+44 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+44 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+44 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+45 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+45 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+45 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+46 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+46 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+46 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+47 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+48 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+48 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+48 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+49 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+49 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+49 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+50 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+50 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+50 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+51 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+51 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+51 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+52 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+52 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+52 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+53 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+53 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+53 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+54 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+54 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+54 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+55 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+55 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+55 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+56 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+56 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+56 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+57 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+58 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+58 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+58 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+59 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+59 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+59 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+60 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+60 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+60 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+61 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+61 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+61 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+62 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+62 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+62 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+63 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+63 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+63 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+64 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+64 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+64 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+65 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+65 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+65 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+66 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+66 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+66 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+67 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+68 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+68 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+68 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+69 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+69 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+69 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+70 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+70 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+70 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+71 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+71 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+71 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+72 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+72 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+72 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+73 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+73 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+73 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+74 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+74 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+74 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+75 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+75 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+75 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+76 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+76 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+76 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+77 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+78 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+78 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+78 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+79 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+79 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+79 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+80 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+80 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+80 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
+81 -16.40331 -21.3038 42.2878 0.1858529 1.537192 1.5 5.54982
+81 1.211788 1.181785 15.23529 0.01610838 0.9175564 1.5 9.313525
+81 5.288862 4.878821 47.73512 0.0664875 3.627845 1.5 6.534336
+82 -11.68283 -14.94316 32.01583 0.4230731 2.729534 1.5 4.960023
+82 0.1583492 0.03833717 1.029285 0.01767288 0.4925226 1.5 12.13987
+82 15.50322 1.761843 21.46617 0.3592727 2.069925 1.5 5.571526
+83 -18.48018 -18.59019 34.70544 0.2684122 2.47351 1.5 5.369024
+83 1.201787 0.8817548 2.871875 0.06907267 3.059562 1.5 9.548344
+83 6.448978 4.448778 12.00845 0.3001395 1.461217 1.5 6.399999
+84 -11.12945 -11.35947 27.01259 0.2786796 2.612895 1.5 5.597309
+84 2.6686 19.77031 11.44234 0.3040265 3.589068 1.5 6.435422
+84 -0.08834217 -0.2183552 0.7805617 0.01205635 1.112686 1.5 12.85493
+85 -12.48958 -18.27016 10.6147 0.5626736 1.937741 1.5 5.849811
+85 0.4017068 0.3617028 0.8806488 0.0293225 1.187549 1.5 11.76021
+85 10.46938 7.349068 11.99595 0.3247444 2.649565 1.5 6.313781
+86 -1.361803 -12.77294 52.24541 0.03493096 3.154616 1.5 7.140002
+86 1.06844 23.03064 44.62869 0.03248871 2.939099 1.5 7.389781
+86 -0.07167383 -0.6217288 2.05649 0.004406628 2.591178 1.5 12.90007
+87 24.45078 35.45188 12.03088 0.3008483 1.463511 1.5 6.390063
+88 -7.239057 -4.318765 17.04554 0.09688595 1.487076 1.5 7.248989
+88 -1.171784 -0.7817448 3.348123 0.02047433 2.902479 1.5 10.70364
+88 18.2835 17.54342 56.52159 0.07901845 2.380442 1.5 6.168811
+89 -12.27289 23.86072 20.47096 0.136267 2.27813 1.5 6.680109
+89 4.452112 -6.238957 28.01408 0.04784984 2.242146 1.5 7.475789
+89 0.8684202 -1.301797 8.005488 0.009684243 5.314029 1.5 10.5703
+90 -0.1216788 0.6784012 1.905763 0.01315577 4.78328 1.5 11.78872
+90 -2.361903 11.99953 14.26937 0.2734991 2.161853 1.5 6.308281
+90 2.418575 -23.46401 12.01995 0.3007466 4.608408 1.5 6.391425
diff --git a/ChiBenchmark/chiBenchmark/timer.h b/ChiBenchmark/chiBenchmark/timer.h
new file mode 100644
index 0000000..781c8f7
--- /dev/null
+++ b/ChiBenchmark/chiBenchmark/timer.h
@@ -0,0 +1,22 @@
+#pragma once
+//
+#include <sys/time.h>
+extern "C"{
+double myseconds()
+{
+ struct timeval tp;
+ struct timezone tzp;
+ //
+ int i = gettimeofday(&tp,&tzp);
+ //
+ return ( (double) tp.tv_sec + (double) tp.tv_usec * 1.e-6 );
+}
+}
+
+
+static inline unsigned long long cycles()
+{
+ unsigned long long u;
+ __asm__ volatile ("rdtscp;shlq $32,%%rdx;orq %%rdx,%%rax;movq %%rax,%0":"=q"(u)::"%rax", "%rdx", "rcx");
+ return u;
+}
diff --git a/ChiBenchmark/lenstoolchi/examples/Jauzac/mod.par~ b/ChiBenchmark/lenstoolchi/examples/Jauzac/mod.par~
index f920c09..d5ed10d 100644
--- a/ChiBenchmark/lenstoolchi/examples/Jauzac/mod.par~
+++ b/ChiBenchmark/lenstoolchi/examples/Jauzac/mod.par~
@@ -1,171 +1,777 @@
runmode
reference 3 109.3982 37.745778
verbose 0
- NOmass 0 400 0.54 mass.fits
+ NOmass 3 400 0.54 mass.fits
inverse 0 0.12 350
chi2 1
end
grille
nlens 96
nlens_opt 4
end
image
NOarcletstat 6 0 wl.cons
NOsigell 0.3
- multfile 3 imawcs.mul
+ multfile 3 theoretical_images_time.txt
wcs 1
- z_m_limit 1 5 1 3.2 5.0 0.1
- z_m_limit 2 34 1 1.0 5.0 0.1
- z_m_limit 4 7 1 1.5 4.0 0.1
- z_m_limit 5 8 1 2.0 3.5 0.1
- z_m_limit 6 18 1 1.8 4.2 0.1
- z_m_limit 8 16 1 2.6 4.5 0.1
- z_m_limit 9 17 1 2.2 3.7 0.1
- z_m_limit 10 2 1 1.5 4.5 0.1
- z_m_limit 11 20 1 0.6 5.0 0.1
- z_m_limit 12 25 1 2.6 4.3 0.1
- z_m_limit 13 29 1 1.4 2.2 0.1
- z_m_limit 14 31 1 1.6 2.0 0.1
- z_m_limit 15 32 1 2.1 3.4 0.1
- z_m_limit 16 33 1 2.6 4.6 0.1
- z_m_limit 17 37 1 1.6 5.0 0.1
- z_m_limit 18 49 1 2.6 5.0 0.1
- z_m_limit 19 50 1 2.6 5.0 0.1
- z_m_limit 20 52 1 2.4 4.0 0.1
- z_m_limit 23 27 1 0.6 5.0 0.1
- z_m_limit 24 36 1 2.0 4.0 0.1
- z_m_limit 25 39 1 2.7 5.0 0.1
- z_m_limit 26 45 1 2.0 3.4 0.1
- z_m_limit 27 55 1 1.4 3.4 0.1
- z_m_limit 28 56 1 2.2 3.6 0.1
- z_m_limit 29 57 1 1.5 2.0 0.1
- z_m_limit 30 58 1 0.6 7.0 0.1
- z_m_limit 31 59 1 2.0 6.0 0.1
- z_m_limit 33 62 1 1.0 4.2 0.1
- z_m_limit 34 63 1 0.7 3.4 0.1
- z_m_limit 35 64 1 1.0 4.5 0.1
- z_m_limit 36 65 1 3.5 6.5 0.1
- z_m_limit 37 66 1 4.0 8.0 0.1
- z_m_limit 38 67 1 2.0 7.0 0.1
- z_m_limit 39 69 1 2.0 7.0 0.1
- z_m_limit 40 70 1 3.0 5.0 0.1
- z_m_limit 41 71 1 1.5 3.5 0.1
- z_m_limit 42 72 1 1.5 4.5 0.1
- z_m_limit 43 73 1 1.0 4.5 0.1
- z_m_limit 44 74 1 1.5 4.5 0.1
- z_m_limit 45 75 1 1.5 5.5 0.1
- z_m_limit 46 76 1 1.5 5.5 0.1
- sigposArcsec 1.8
forme -1
end
potentiel 1
- profil 81
- x_centre 0.0
- y_centre 0.0
- ellipticite 0.15
- angle_pos 60.0
- core_radius 10.0
- cut_radius 500.0
- v_disp 800.0
- z_lens 0.54
- end
-limit 1
- x_centre 1 -10. 10. 1.0
- y_centre 1 -15. 0. 1.0
- ellipticite 1 0.05 0.4 0.01
- angle_pos 1 40.0 90. 1.0
- core_radius 1 8.0 25.0 1.
- cut_radius 1 100.0 1000.0 1.
- v_disp 1 750. 950. 1.
- end
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 800. #Dispersion Velocity [km/s]
+ rcut 500 #Cut radius (PIEMD distribution)
+ rcore 5 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
potentiel 2
- profil 81
- x_centre 35.0
- y_centre -10.0
- ellipticite 0.15
- angle_pos 60.0
- core_radius 5.0
- cut_radius 500.0
- v_disp 600.0
- z_lens 0.54
- end
-limit 2
- x_centre 1 20. 30. 1.0
- y_centre 1 -29. -15. 1.0
- ellipticite 1 0.3 0.7 0.01
- angle_pos 1 35.0 70. 1.0
- core_radius 1 1.0 12.0 1.
- cut_radius 1 100.0 1000.0 1.
- v_disp 1 350. 600. 1.
- end
-potentiel 3
- profil 81
- x_centre 65.0
- y_centre 40.0
- ellipticite 0.15
- angle_pos 60.0
- core_radius 10.0
- cut_radius 500.0
- v_disp 700.0
- z_lens 0.54
- end
-limit 3
- x_centre 1 30. 85. 1.0
- y_centre 1 30. 55. 1.0
- ellipticite 1 0.3 0.6 0.01
- angle_pos 1 0.0 180. 1.0
- core_radius 1 10.0 30.0 1.
- cut_radius 1 100.0 1000.0 1.
- v_disp 1 650. 950. 1.
- end
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 25.000 #X Position [arcsec]
+ y_centre 20.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 140. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
potentiel 4
- profil 81
- x_centre 130.0
- y_centre 80.0
- ellipticite 0.15
- angle_pos 60.0
- core_radius 10.0
- cut_radius 500.0
- v_disp 1000.0
- z_lens 0.54
- end
-limit 4
- x_centre 1 110. 135. 1.0
- y_centre 1 60. 80. 1.0
- ellipticite 1 0.45 0.7 0.01
- angle_pos 1 0.0 180. 1.0
- core_radius 1 10.0 30.0 1.
- cut_radius 1 100.0 1000.0 1.
- v_disp 1 800. 1000. 1.
- end
-potentiel 683
- profil 81
- x_centre -19.4
- y_centre -21.7
- ellipticite 0.23
- angle_pos -40.0
- core_radius 1.0
- cut_radius 150.0
- v_disp 180.0
- z_lens 0.15
- end
-potfile galaxies
- filein 1 gal2.cat
- zlens 0.54
- profil 81
- corekpc 0.3
- mag0 20.66
- sigma 1 150. 350.
- cutkpc 1 45. 70.
- slope 0 4.
- vdslope 0 4.
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre -16.000 #X Position [arcsec]
+ y_centre 2.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 130. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 5
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre -10.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 110. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 6
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 3.000 #X Position [arcsec]
+ y_centre -10.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 80. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 7
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 17.000 #X Position [arcsec]
+ y_centre -7.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 120. #Dispersion Velocity [km/s]
+ rcut 150 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 8
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 9
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 10
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 11
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 12
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 13
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 14
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 15
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 16
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 17
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 18
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 19
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 20
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 21
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 22
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 23
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 24
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 25
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 26
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 27
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 28
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 29
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 30
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 8
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 9
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 10
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 11
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 12
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 13
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 14
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 15
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 16
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 17
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 18
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 19
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 20
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 21
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 22
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 23
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 24
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 25
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 26
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 27
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 28
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 29
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 30
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 8
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 9
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 10
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 11
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 12
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 13
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 14
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 15
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 16
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 17
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 18
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 19
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 20
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 21
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 22
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 23
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 24
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 25
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 26
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 27
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 28
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 29
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
+ end
+potentiel 30
+ profil 8 #Profile: 5 SIS, 8 PIEMD
+ x_centre 0.000 #X Position [arcsec]
+ y_centre 0.000 #Y Position [arcsec]
+ ellipticity 0.11 #Ellipticity
+ v_disp 1. #Dispersion Velocity [km/s]
+ rcut 5 #Cut radius (PIEMD distribution)
+ rcore 1 #Core radius (PIEMD distribution)
+ z_lens 0.4 #Redshift of lens
end
cosmology
H0 70.0
omega 0.3
lambda 0.7
end
champ
dmax 400.
end
fini