Commit c98c6560 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden
Browse files

Ring find in extrinsic burning

parent bd2b0518
Pipeline #163099 passed with stages
in 2 minutes and 17 seconds
......@@ -35,7 +35,6 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include <lal/TimeDelay.h>
#include <lal/DetectorSite.h>
#include <lal/DetResponse.h>
// #include "BayesWave.h"
/*******************************************************************************************
......@@ -216,8 +215,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
alpha = gsl_rng_uniform(r);
// if(alpha > 0.80 && net->Nifo > 1) // ring
if(alpha > 1.0 && net->Nifo > 1) // ring
if(alpha > 0.8 && net->Nifo > 1) // ring
{
......@@ -240,7 +238,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
// get new sky locations with the same time delays
Ring(net, skyx[q], skyy[q], id1, id2, r, rundata->gmst);
// for(i = 0; i < NS; i++) skyy[q][i] = skyx[q][i];
qyx = 1.0;
qxy = 1.0;
......@@ -357,42 +355,42 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
} // ends choice of update
/*
if(mc%100 == 0)
{
ic = who[1];
phi = skyx[ic][0];
if(phi > TPI) phi -= TPI;
if(phi < 0.0) phi += TPI;
skyx[ic][0] = phi;
Mchirp = exp(paramx[ic][0])/MSUN_SI;
Mtot = exp(paramx[ic][1])/MSUN_SI;
eta = pow((Mchirp/Mtot), (5.0/3.0));
dm = sqrt(1.0-4.0*eta);
m1 = Mtot*(1.0+dm)/2.0;
m2 = Mtot*(1.0-dm)/2.0;
chieff = (m1*paramx[ic][2]+m2*paramx[ic][3])/Mtot;
// if(mc%100 == 0)
// {
// ic = who[1];
// phi = skyx[ic][0];
// if(phi > TPI) phi -= TPI;
// if(phi < 0.0) phi += TPI;
// skyx[ic][0] = phi;
// Mchirp = exp(paramx[ic][0])/MSUN_SI;
// Mtot = exp(paramx[ic][1])/MSUN_SI;
// eta = pow((Mchirp/Mtot), (5.0/3.0));
// dm = sqrt(1.0-4.0*eta);
// m1 = Mtot*(1.0+dm)/2.0;
// m2 = Mtot*(1.0-dm)/2.0;
// chieff = (m1*paramx[ic][2]+m2*paramx[ic][3])/Mtot;
// counter, log likelihood, chirp mass, total mass, effective spin, phase shift , time shift, distance, RA , sine of DEC,
// polarization angle, cos inclination
// // counter, log likelihood, chirp mass, total mass, effective spin, phase shift , time shift, distance, RA , sine of DEC,
// // polarization angle, cos inclination
DL = exp(pallx[ic][6])/(1.0e6*PC_SI*skyx[ic][4]);
z = z_DL(DL);
// Note that skyx[ic][5], skyx[ic][6], hold different quantities than what is printed by the other MCMCs
// DL = exp(pallx[ic][6])/(1.0e6*PC_SI*skyx[ic][4]);
// z = z_DL(DL);
// // Note that skyx[ic][5], skyx[ic][6], hold different quantities than what is printed by the other MCMCs
fprintf(chain,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", mxc[1], logLx[ic], Mchirp, Mtot, chieff, skyx[ic][5], \
skyx[ic][6], DL, skyx[ic][0], skyx[ic][1], \
skyx[ic][2], skyx[ic][3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
// // fprintf(rundata->chainS,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\n", mxc[1], logLx[ic], Mchirp, Mtot, chieff, skyx[ic][5], \
// skyx[ic][6], DL, skyx[ic][0], skyx[ic][1], \
// skyx[ic][2], skyx[ic][3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
mxc[1] += 1;
// mxc[1] += 1;
} */
// }
if(mc%10000 == 0)
{
......@@ -1030,16 +1028,39 @@ void skymap(struct Net *net, double *paramsx, double *paramsy, int ifo1, int ifo
// extract the primitives f+, fx for each detector at x
alpha = paramsx[0];
sindelta = paramsx[1];
ComputeDetFant(net, 0.0, alpha, sindelta, &fxplus[1], &fxcross[1], ifo1, gmst);
ComputeDetFant(net, 0.0, alpha, sindelta, &fxplus[2], &fxcross[2], ifo2, gmst);
// Detector 1
ComputeDetFant(net, 0.0, alpha, sindelta, &x, &y, ifo1, gmst);
psi0 = 0.5*atan2(y,x);
ComputeDetFant(net, psi0, alpha, sindelta, &x, &y, ifo1, gmst);
fxplus[1] = x*cos(2.0*psi0);
fxcross[1] = fxplus[1]*tan(2.0*psi0);
// Detector 2
ComputeDetFant(net, 0.0, alpha, sindelta, &x, &y, ifo2, gmst);
psi0 = 0.5*atan2(y,x);
ComputeDetFant(net, psi0, alpha, sindelta, &x, &y, ifo2, gmst);
fxplus[2] = x*cos(2.0*psi0);
fxcross[2] = fxplus[2]*tan(2.0*psi0);
TimeDelays(net, alpha, sindelta, dtimesx, gmst);
// extract the primitives f+, fx for each detector at y
alpha = paramsy[0];
sindelta = paramsy[1];
ComputeDetFant(net, 0.0, alpha, sindelta, &fyplus[1], &fycross[1], ifo1, gmst);
ComputeDetFant(net, 0.0, alpha, sindelta, &fyplus[2], &fycross[2], ifo2, gmst);
// Detector 1
ComputeDetFant(net, 0.0, alpha, sindelta, &x, &y, ifo1, gmst);
psi0 = 0.5*atan2(y,x);
ComputeDetFant(net, psi0, alpha, sindelta, &x, &y, ifo1, gmst);
fyplus[1] = x*cos(2.0*psi0);
fycross[1] = fyplus[1]*tan(2.0*psi0);
// Detector 2
ComputeDetFant(net, 0.0, alpha, sindelta, &x, &y, ifo2, gmst);
psi0 = 0.5*atan2(y,x);
ComputeDetFant(net, psi0, alpha, sindelta, &x, &y, ifo2, gmst);
fyplus[2] = x*cos(2.0*psi0);
fycross[2] = fyplus[2]*tan(2.0*psi0);
TimeDelays(net, alpha, sindelta, dtimesy, gmst);
......@@ -1068,6 +1089,7 @@ void skymap(struct Net *net, double *paramsx, double *paramsy, int ifo1, int ifo
}
void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng * r, double gmst)
{
......@@ -1082,7 +1104,6 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
double calpha, salpha;
double H1mag, L1mag, zmag;
double cmu, smu;
double gmstrad;
H1 = (double*)malloc(sizeof(double)* 3);
L1 = (double*)malloc(sizeof(double)* 3);
......@@ -1161,18 +1182,16 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
//remap gmst back to [0:2pi]
double intpart;
double decpart;
gmstrad = gmst/TWOPI;
double gmstrad = gmst/TWOPI;
intpart = (int)( gmstrad );
decpart = gmstrad - (double)intpart;
gmstrad = decpart*TWOPI;
double dec;
alpha = skyx[0];
sind = skyx[1];
cosd = sqrt(1.0-sind*sind);
dec = asin(sind);
calpha = cos(gmstrad-alpha);
salpha = -sin(gmstrad-alpha);
calpha = cos(alpha-gmstrad);
salpha = -sin(alpha-gmstrad);
// Unusual convention on spherical angle
// propagation direction
......@@ -1180,6 +1199,7 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
kv[1] = cosd*salpha;
kv[2] = sind;
// cosine of angle between propagation direction and vector connecting the sites
cmu = 0.0;
for(i = 0; i <3; i++) cmu += zv[i]*kv[i];
......@@ -1202,21 +1222,15 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
rot = TPI*gsl_rng_uniform(r);
cr = cos(rot);
sr = sin(rot);
double c1momega = 1.0 - rot;
// find rotated k vector
// for(i = 0; i <3; i++) kvr[i] = cmu*zv[i] + smu*(cr*uv[i]+sr*vv[i]);
kvr[0] = (c1momega*zv[0]*zv[0] + cr) *kv[0] + (c1momega*zv[0]*zv[1] - sr*zv[2])*kv[1] + (c1momega*zv[0]*zv[2] + sr*zv[1])*kv[2];
kvr[1] = (c1momega*zv[0]*zv[1] + sr*zv[2])*kv[0] + (c1momega*zv[1]*zv[1] + cr) *kv[1] + (c1momega*zv[1]*zv[2] - sr*zv[0])*kv[2];
kvr[2] = (c1momega*zv[0]*zv[2] - sr*zv[1])*kv[0] + (c1momega*zv[1]*zv[2] + sr*zv[0])*kv[1] + (c1momega*zv[2]*zv[2] + cr) *kv[2];
for(i = 0; i <3; i++) kvr[i] = cmu*zv[i] + smu*(cr*uv[i]+sr*vv[i]);
skyy[1] = kvr[2];
rot = atan2(kvr[1],kvr[0]);
rot = gmstrad + rot;
rot = gmstrad - rot;
if(rot < 0.0) rot += TPI;
if(rot > TPI) rot -= TPI;
skyy[0] = rot;
......@@ -1231,154 +1245,6 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
}
// void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng * r, double gmst)
// {
// int i;
// double *H1;
// double *L1;
// double *V1;
// double *kv, *kvr;
// double *zv, *uv, *vv;
// double sind, cosd, rot, alpha, x;
// double cr, sr;
// double calpha, salpha;
// double H1mag, L1mag, zmag;
// double cmu, smu;
// H1 = (double*)malloc(sizeof(double)* 3);
// L1 = (double*)malloc(sizeof(double)* 3);
// V1 = (double*)malloc(sizeof(double)* 3);
// zv = (double*)malloc(sizeof(double)* 3);
// uv = (double*)malloc(sizeof(double)* 3);
// vv = (double*)malloc(sizeof(double)* 3);
// kv = (double*)malloc(sizeof(double)* 3);
// kvr = (double*)malloc(sizeof(double)* 3);
// // From https://dcc.ligo.org/public/0072/P000006/000/P000006-C.pdf (note sign flip on y axis)
// // Location of Hanford vertex from Geocenter in s
// H1[0] = net->location[0][0]; // -2.161414928e6/CLIGHT;
// H1[1] = net->location[0][1]; //3.834695183e6/CLIGHT;
// H1[2] = net->location[0][2]; // 4.600350224e6/CLIGHT;
// //H1mag = sqrt(H1[0]*H1[0]+H1[1]*H1[1]+H1[2]*H1[2]);
// // Location of Livingston vertex from Geocenter in s
// L1[0] = net->location[1][0]; // -74276.04192/CLIGHT;
// L1[1] = net->location[1][1]; //5.496283721e6/CLIGHT;
// L1[2] = net->location[1][2]; // 3.224257016e6/CLIGHT;
// //L1mag = sqrt(L1[0]*L1[0]+L1[1]*L1[1]+L1[2]*L1[2]);
// // Location of Virgo vertex from Geocenter in s
// V1[0] = net->location[2][0]; //4546374.098/CLIGHT;
// V1[1] = net->location[2][1]; //-842989.6972/CLIGHT;
// V1[2] = net->location[2][2]; // 4378576.963/CLIGHT;
// if(d1 == 0)
// {
// if(d2 == 1)
// {
// for(i = 0; i <3; i++) zv[i] = H1[i]-L1[i];
// }
// if(d2 == 2)
// {
// for(i = 0; i <3; i++) zv[i] = H1[i]-V1[i];
// }
// }
// if(d1 == 1)
// {
// if(d2 == 0)
// {
// for(i = 0; i <3; i++) zv[i] = L1[i]-H1[i];
// }
// if(d2 == 2)
// {
// for(i = 0; i <3; i++) zv[i] = L1[i]-V1[i];
// }
// }
// if(d1 == 2)
// {
// if(d2 == 0)
// {
// for(i = 0; i <3; i++) zv[i] = V1[i]-H1[i];
// }
// if(d2 == 1)
// {
// for(i = 0; i <3; i++) zv[i] = V1[i]-L1[i];
// }
// }
// zmag = 0.0;
// for(i = 0; i <3; i++) zmag += zv[i]*zv[i];
// zmag = sqrt(zmag);
// // Unet->Nifot vector between the sites
// for(i = 0; i <3; i++) zv[i] /= zmag;
// alpha = skyx[0];
// sind = skyx[1];
// cosd = sqrt(1.0-sind*sind);
// calpha = cos(alpha-gmst);
// salpha = -sin(alpha-gmst);
// // Unusual convention on spherical angle
// // propagation direction
// kv[0] = cosd*calpha;
// kv[1] = cosd*salpha;
// kv[2] = sind;
// // cosine of angle between propagation direction and vector connecting the sites
// cmu = 0.0;
// for(i = 0; i <3; i++) cmu += zv[i]*kv[i];
// if(cmu > 0.999999) cmu = 0.999999;
// if(cmu < -0.999999) cmu = -0.999999;
// smu = sqrt(1.0-cmu*cmu);
// // vector in plane of prop vector and line connecting sites
// for(i = 0; i <3; i++) uv[i] = (kv[i] - cmu*zv[i])/smu;
// //completing the triad
// vv[0] = (zv[1]*kv[2]-zv[2]*kv[1])/smu;
// vv[1] = (zv[2]*kv[0]-zv[0]*kv[2])/smu;
// vv[2] = (zv[0]*kv[1]-zv[1]*kv[0])/smu;
// // Rotating about zv axis will produce new sky locations with the same time delay between H1 and L1
// rot = TPI*gsl_rng_uniform(r);
// cr = cos(rot);
// sr = sin(rot);
// // find rotated k vector
// for(i = 0; i <3; i++) kvr[i] = cmu*zv[i] + smu*(cr*uv[i]+sr*vv[i]);
// skyy[1] = kvr[2];
// rot = atan2(kvr[1],kvr[0]);
// rot = gmst - rot;
// if(rot < 0.0) rot += TPI;
// skyy[0] = rot;
// free(H1);
// free(L1);
// free(V1);
// free(kv);
// free(kvr);
// free(zv);
// free(uv);
// free(vv);
// }
/* ********************************************************************************** */
/* */
/* Wrappers for LAL */
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment