#Problem1

input=scan(file = "sn_data_riess.dat", what = list(character(), double(), double(), double()), skip=1, multi.line=FALSE)
data=cbind(input[[2]], input[[3]], input[[4]])

#Hubble expansion in Dark Energy models
#we assume flat Universe
H0=72.0 #km/s/Mpc
c=3*10^5 #km/s
N=186 #number of data

#################Model_A########################
#Hubble law for LCDM
Hubble_A=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))

#H^(-1)
Hinv_A=function(zint,Omegam) 1.0/Hubble_A(zint,Omegam)

#Luminosity distance
dL_A=function(z,Omegami){ 
	dLsol=z
	for (i in (1:length(z))){
		zarg=z[i]
		I=integrate(Hinv_A,0.0,zarg,Omegam=Omegami)
		dLsol[i] = c*(1+zarg)*I$value
		}
	return(dLsol)
	}

#Magnitude
mag_A=function(z,Omegam,M) M+5.0*log10(H0*dL_A(z,Omegam))

# Define function chi2 
chi2_A=function(Omegam,M){
chi2sol=0.0
for (i in (1:N))  
               {          
		chi2sol=chi2sol+((data[i,2]-mag_A(data[i,1],Omegam,M))^2)/(data[i,3])^2
	       }
	return(chi2sol)
}

###############Model_B####################
#Hubble law for wCDM
Hubble_B=function(z,Omegam,w) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam)*(1+z)^(3*(1+w)))

#H^(-1)
Hinv_B=function(zint,Omegam, w) 1.0/Hubble_B(zint,Omegam, w)


#Luminosity distance
dL_B=function(z,Omegami, wi){ 
	dLsol=z
	for (i in (1:length(z))){
		zarg=z[i]
		I=integrate(Hinv_B,0.0,zarg,Omegam=Omegami, w=wi)
		dLsol[i] = c*(1+zarg)*I$value
		}
	return(dLsol)
	}

#Magnitude
mag_B=function(z,Omegam,M, w) M+5.0*log10(H0*dL_B(z,Omegam, w))

# Define function chi2 
chi2_B=function(Omegam,M ,w){
chi2sol=0.0
for (i in (1:N))  
               {          
		chi2sol=chi2sol+((data[i,2]-mag_B(data[i,1],Omegam,M,w))^2)/(data[i,3])^2
	       }
	return(chi2sol)
}


#i) Omegam=0,3, w=m{-2.0, -1.5, -1.0, -0.5}
plot(data[,1], dL_B(data[,1], 0.3, -2.0),type="p", col="black", pch=19, cex=.6, xlab='redshift z',ylab='lumiosity distance dL(z)')
points(data[,1], dL_B(data[,1], 0.3, -1.5),type="p", col="blue", pch=19, cex=.6)
points(data[,1], dL_B(data[,1], 0.3, -1.0),type="p", col="green", pch=19, cex=.6)
points(data[,1], dL_B(data[,1], 0.3, -0.5),type="p", col="red", pch=19, cex=.6)

X11()
#ii) w=-1, Omegam={0.1, 0.2, 0.3, 0.4, 0.5}
plot(data[,1], dL_B(data[,1], 0.1, -1.0),type="p", col="black", pch=19, cex=.6, xlab='redshift z',ylab='lumiosity distance dL(z)')
points(data[,1], dL_B(data[,1], 0.2, -1.0),type="p", col="blue", pch=19, cex=.6)
points(data[,1], dL_B(data[,1], 0.3, -1.0),type="p", col="green", pch=19, cex=.6)
points(data[,1], dL_B(data[,1], 0.4, -1.0),type="p", col="red", pch=19, cex=.6)
points(data[,1], dL_B(data[,1], 0.5, -1.0),type="p", col="yellow", pch=19, cex=.6)


#Problem 2


X11()
#i) Model_A: Omegam variiert von 0.1 bis 0.5
plot(data[,1], mag_A(data[,1], 0.1, 16),type="p", col="black", pch=19, cex=.6, xlab='redshift z',ylab='theoretical magnitude m')
points(data[,1], mag_A(data[,1], 0.2, 16),type="p", col="blue", pch=19, cex=.6)
points(data[,1], mag_A(data[,1], 0.3, 16),type="p", col="green", pch=19, cex=.6)
points(data[,1], mag_A(data[,1], 0.4, 16),type="p", col="red", pch=19, cex=.6)
points(data[,1], mag_A(data[,1], 0.5, 16),type="p", col="yellow", pch=19, cex=.6)

#ii)
X11()
# Model_B: Omegam variiert von 0.1 bis 0.5 bei w=-1.0
plot(data[,1], mag_B(data[,1], 0.1, 16, -1.0),type="p", col="black", pch=19, cex=.6, xlab='redshift z',ylab='theoretical magnitude m')
points(data[,1], mag_B(data[,1], 0.2, 16, -1.0),type="p", col="blue", pch=19, cex=.6)
points(data[,1], mag_B(data[,1], 0.3, 16, -1.0),type="p", col="green", pch=19, cex=.6)
points(data[,1], mag_B(data[,1], 0.4, 16, -1.0),type="p", col="red", pch=19, cex=.6)
points(data[,1], mag_B(data[,1], 0.5, 16, -1.0),type="p", col="yellow", pch=19, cex=.6)

X11()
# Model_B: w variiert von -2.0 bis -0.5 bei Omegam=0.3
plot(data[,1], mag_B(data[,1], 0.3, 16, -2.0),type="p", col="black", pch=19, cex=.6, xlab='redshift z',ylab='theoretical magnitude m')
points(data[,1], mag_B(data[,1], 0.3, 16, -1.5),type="p", col="blue", pch=19, cex=.6)
points(data[,1], mag_B(data[,1], 0.3, 16, -1.0),type="p", col="green", pch=19, cex=.6)
points(data[,1], mag_B(data[,1], 0.3, 16, -0.5),type="p", col="red", pch=19, cex=.6)

#Problem 3
#Model_A: chi2 for Omegam={0.0,1.0} and M={15.5, 16.5}
values_A=10
Omegam_vec_A=seq(0.0, 1.0, by = 1/values_A)
M_vec_A=seq(15.5, 16.5, by=1/values_A)
twodim=array(0, dim=c(values_A+1,values_A+1))
for (i in (0:values_A+1)){
	for (j in (0:values_A+1)){
		twodim[i,j]=chi2_A(Omegam_vec_A[i], M_vec_A[j])
	}
}

#Model_B: chi2 for Omegam={0.0,1.0} and M={15.5, 16.5}, w={-2.0,-0.5}
values_B=10
Omegam_vec_B=seq(0.0, 1.0, by = 1/values_B)
M_vec_B=seq(15.5, 16.5, by=1/values_B)
w_vec=seq(-2.0, -0.5, by=1/(values_B/1.5))

threedim=array(0, dim=c(values_B+1,values_B+1, values_B+1))
for (i in (0:values_B+1)){
	for (j in (0:values_B+1)){
		for (k in (0:values_B+1)){
			threedim[i,j,k]=chi2_B(Omegam_vec_B[i], M_vec_B[j], w_vec[k])
		}
	}
}
#Find best fit parameters: where is the minimum of chi2?
#Model_A: find minimum
min_A=twodim[1,1]
min_A_pos=array(0, dim=c(2))
for (i in (0:values_A+1)){
	for (j in (0:values_A+1)){
		if (twodim[i,j] < min_A){
			min_A=twodim[i,j]
			min_A_pos[1]=i
			min_A_pos[2]=j
		}
	}
}
twodim
out1=c("chi2min",min_A)
out2=c("Omegam_min",Omegam_vec_A[min_A_pos[1]])
out3=c("M_min",M_vec_A[min_A_pos[2]])                                
print(out1)
print(out2)
print(out3)

#Model_B: find minimum
min_B=threedim[1,1,1]
min_B_pos=array(0, dim=c(3))
for (i in (0:values_B+1)){
	for (j in (0:values_B+1)){
		for (k in (0:values_B+1)){
		
			if (threedim[i,j,k] < min_B){
				min_B=threedim[i,j,k]
				min_B_pos[1]=i
				min_B_pos[2]=j
				min_B_pos[3]=k
			}
		}
	}
}
threedim
out1=c("chi2min",min_B)
out2=c("Omegam_min",Omegam_vec_B[min_B_pos[1]])
out3=c("M_min",M_vec_B[min_B_pos[2]]) 
out4=c("w_min",w_vec[min_B_pos[3]])                               
print(out1)
print(out2)
print(out3)
print(out4)

#Problem4

#Model_A posterior function  
posterior_A=array(0, c(values_A+1, values_A+1)) 
for (i in (0:values_A)){
	for (j in (0:values_A)){
		posterior_A[i,j]=exp(-0.5*(twodim[i,j]-min_A))
	}
}                    
  
#Model_B posterior function                        
posterior_B=array(0, c(values_B+1, values_B+1, values_B+1))
for (i in (0:values_B)){
	for (j in (0:values_B)){
		for (k in (0:values_B)){
			posterior_B[i,j,k]=exp(-0.5*(threedim[i,j,k]-min_B))
		}
	}
}    

give_function=function(i,j) f=splinefun(w_vec, posterior_B[i,j,])
Myintegrate = function(F,LOW,UP){
    pI = integrate(F,LOW,UP)
    return(pI$value)
}

integration_B=array(0, c(values_B+1, values_B+1))
for (i in (1:values_B)){
	for (j in (1:values_B)){
		I=Myintegrate(give_function(i,j),-2.0,-0.5)
		integration_B[i,j]=I
	}
}

X11()
contour(Omegam_vec_A,M_vec_A,posterior_A,drawlabels=FALSE,xlab='Omegam',ylab='M',xlim=c(0,1),ylim=c(15.5,16.5),col = "red")
contour(Omegam_vec_B,M_vec_B,integration_B,add = TRUE)

#Problem 5
# Marginalize over M


give_function_MA=function(i) f=splinefun(M_vec_A, posterior_A[i,])
give_function_MB=function(i,ii) f=splinefun(M_vec_B, posterior_B[i,,ii])


integration_MA=array(0, c(values_A+1))
for (i in (1:values_A)){
		I=Myintegrate(give_function_MA(i),15.5,16.5)
		integration_MA[i]=I
	}

integration_MB=array(0, c(values_B+1, values_B+1))
for (i in (1:values_B)){
	for (ii in (1:values_B)){
		I=Myintegrate(give_function_MB(i,ii),15.5,16.5)
		integration_MB[i,ii]=I
	}
}
X11()
contour(Omegam_vec_B,w_vec,integration_MB,xlab='Omegam',ylab='w',col = "blue")

#The graphical result of problem 4 shows, that with wCDM theory you get a smaller error for the same resolution (values_A=value_B=10).
#This is supported by the result of problem 3, where the chi2min of LambdaCDM is 238, the chi2min of wCDM is 229, so that the latter
#better fits to the data.