QdotSupernova[tt_]:=(
c=3*(10^10);
Msun=2*(10^33);
λN=(0.693/(6.075*24*60*60));
λC=(0.693/(77.236*24*60*60));
EN=1.72*(10^6)*1.6*(10^(-19))*(10^7);
EC=3.57*(10^6)*1.6*(10^(-19))*(10^7);
N0=0.01*(Msun/((59)))*6*(10^23);
E0=(10^51);



Qdot[tt]=λN*N0(((λC*EC)/(λC-λN))*(Exp[-λN*tt]-Exp[-λC*tt])+EN*Exp[-λN*tt]))

Lum[t_,Mej_,Vej_,kappa_,E0_,Rin_]:=(c=3*(10^10);κ=kappa;v=Vej*c;
Msun=2*(10^33);
M=Mej*Msun;
Rsun=6.955*(10^10);


Tdiff=(Sqrt[3κ*(M/(4π c*v))]);
R0=Rin;
Rdiff=v*Tdiff;
InEn[t]:=(R0/Rdiff)*(E0/Tdiff)*Exp[((-(t^2))/(2*(Tdiff^2)))];

Qdot[t]=QdotSupernova[t];


Qdotkilonava=0.5*(10^10)*(((t/86400))^(-1.3))*M;

SecEn[t]:=((Exp[((-(t^2))/(2*(Tdiff^2)))])/(Tdiff^2))*Integrate[T*Exp[((T^2)/(2*(Tdiff^2)))]*(*Qdot[t]*)Qdotkilonava,{T,0,t}];
Lum[t]=InEn[t]+SecEn[t])

qn. 3

E0=0, No R0 is needed

Plot[Lum[t,3,0.01,0.1,0,0],{t,200,50000000},AxesLabel->{"time(s)","Luminosity(erg/s)"}]

Qn. 4

E0=(10^51)

LogPlot[{Lum[t,3,0.01,0.1,0,0],Lum[t,3,0.01,0.1,(10^51),Rsun]},{t,200,50000000},AxesLabel->{"time(s)","Luminosity(erg/s)"},PlotLegends->{"qn3","qn4"}]

Qn5

Teff[mm_,vej_,kappa_]:=(c=3*(10^10);κ=kappa;v=vej*c;
Msun=2*(10^33);
M=mm*Msun;Lpeakvalue=(10^40)*((((M/(0.01*Msun))))^(1/3))*((((v/(0.2*c))))^(2/3))*((((κ/10)))^(-2/3));

Tdiff=(Sqrt[3κ*(M/(4π c*v))]);
R0=Rin;
Rdiff=v*Tdiff;
σ=5.5*(10^(-5));

lpeak=(10^40)*((((M/(0.01*Msun))))^(1/3))*((((v/(0.02*c))))^(2/3))*((((κ/10)))^(-2/3));


teff=(((((Lum[t,3,0.01,0.1,0,0])/(4π σ*(Rdiff^2)))))^(1/4)))

Plot[Teff[t,0.3,0.1],{t,200,(10^7)},AxesLabel->{"Time(s)","Eff-Tem(K)"}]

kilonova=Import["/Acdmc/GWs/GWs_summerSchool_2018/Kenta/obs_kilonova.dat"];

LogPlot[Lum[t,0.03,0.2,0.1,(10^51),(10^7)],{t,(10^4),(10^7)}]

(*the 9th column where the t-t0 and also in days. Multiply with 86400 to get it in seconds*)

bnslum=Table[{kilonova[[ii,9]]*86400,kilonova[[ii,2]]},{ii,2,Length[kilonova]}];

Lum[(10^5),0.03,0.2,0.1,(10^51),(10^7)]