diff --git a/IceRayTracing.cc b/IceRayTracing.cc new file mode 100644 index 0000000..98ab073 --- /dev/null +++ b/IceRayTracing.cc @@ -0,0 +1,3039 @@ +/* + This is the IceRayTracing namespace. Author: Uzair Latif + released under GPL3. +*/ + +#include "IceRayTracing.hh" + +void IceRayTracing::SetA(double &A){ + IceRayTracing::A_ice=A; +} + +void IceRayTracing::SetB(double &B){ + IceRayTracing::B_ice=B; +} + +void IceRayTracing::SetC(double &C){ + IceRayTracing::C_ice=C; +} + +/* Get the value of the B parameter for the refractive index model */ +double IceRayTracing::GetB(double z){ + z=fabs(z); + double B=0; + + B=IceRayTracing::B_ice; + + // if(z<=IceRayTracing::TransitionBoundary){ + // B=-0.5019; + // }else{ + // B=-0.448023; + // } + + return B; +} + +/* Get the value of the C parameter for the refractive index model */ +double IceRayTracing::GetC(double z){ + z=fabs(z); + double C=0; + + C=IceRayTracing::C_ice; + + // if(z<=IceRayTracing::TransitionBoundary){ + // C=0.03247; + // }else{ + // C=0.02469; + // } + + return C; +} + + +/* Get the value of refractive index model for a given depth */ +double IceRayTracing::Getnz(double z){ + z=fabs(z); + return IceRayTracing::A_ice+IceRayTracing::GetB(z)*exp(-IceRayTracing::GetC(z)*z); +} + +/* E-feild Power Fresnel coefficient for S-polarised wave which is perpendicular to the plane of propogation/incidence. This function gives you back the reflectance. The transmittance is T=1-R */ +double IceRayTracing::Refl_S(double thetai){ + + double Nair=1; + double Nice=IceRayTracing::Getnz(0); + double n1=Nice; + double n2=Nair; + + double sqterm=sqrt(1-pow((n1/n2)*(sin(thetai)),2)); + double num=n1*cos(thetai)-n2*sqterm; + double den=n1*cos(thetai)+n2*sqterm; + double RS=(num*num)/(den*den); + + if(std::isnan(RS)){ + RS=1; + } + return (RS); +} + +/* E-feild Power Fresnel coefficient for P-polarised wave which is parallel to the plane of propogation/incidence. This function gives you back the reflectance. The transmittance is T=1-R */ +double IceRayTracing::Refl_P(double thetai){ + + double Nair=1; + double Nice=IceRayTracing::Getnz(0); + double n1=Nice; + double n2=Nair; + + double sqterm=sqrt(1-pow((n1/n2)*(sin(thetai)),2)); + double num=n1*sqterm-n2*cos(thetai); + double den=n1*sqterm+n2*cos(thetai); + double RP=(num*num)/(den*den); + if(std::isnan(RP)){ + RP=1; + } + return (RP); +} + +/* The temperature and attenuation model has been taken from AraSim which also took it from here http://icecube.wisc.edu/~araproject/radio/ . This is basically Matt Newcomb's icecube directory which has alot of information, plots and codes about South Pole Ice activities. Please read it if you find it interesting. */ + +/* Temperature model:The model takes in value of depth z in m and returns the value of temperature in Celsius.*/ +double IceRayTracing::GetIceTemperature(double z){ + double depth=fabs(z); + double t = 1.83415e-09*pow(depth,3) + (-1.59061e-08*pow(depth,2)) + 0.00267687*depth + (-51.0696 ); + return t; +} + +/* Ice Attenuation Length model: Takes in value of frequency in Ghz and depth z and returns you the value of attenuation length in m */ +double IceRayTracing::GetIceAttenuationLength(double z, double frequency){ + + double t =IceRayTracing::GetIceTemperature(z); + const double f0=0.0001, f2=3.16; + const double w0=log(f0), w1=0.0, w2=log(f2), w=log(frequency); + const double b0=-6.74890+t*(0.026709-t*0.000884); + const double b1=-6.22121-t*(0.070927+t*0.001773); + const double b2=-4.09468-t*(0.002213+t*0.000332); + double a,bb; + if(frequency<1.){ + a=(b1*w0-b0*w1)/(w0-w1); + bb=(b1-b0)/(w1-w0); + } + else{ + a=(b2*w1-b1*w2)/(w1-w2); + bb=(b2-b1)/(w2-w1); + } + double Lval=1./exp(a+bb*w); + return Lval; +} + +/* Setup the integrand to calculate the attenuation */ +double IceRayTracing::AttenuationIntegrand (double x, void * params) { + + double *p=(double*)params; + + double A0=p[0]; + double Frequency=p[1]; + double L=p[2]; + + double Integrand=(A0/IceRayTracing::GetIceAttenuationLength(x,Frequency))*sqrt(1+pow(tan(asin(L/IceRayTracing::Getnz(x))) ,2)); + return Integrand; +} + +/* Integrate over the integrand to calculate the attenuation */ +double IceRayTracing::IntegrateOverLAttn (double A0, double Frequency, double z0, double z1, double Lvalue) { + gsl_integration_workspace * w= gsl_integration_workspace_alloc (1000); + + double result, error; + double zlimit[2] = {z0,z1}; + double param[3] = {A0,Frequency,Lvalue}; + + gsl_function F; + F.function = &AttenuationIntegrand; + F.params = ¶m; + + gsl_integration_qags (&F, zlimit[0], zlimit[1], 0, 1e-7, 1000, + w, &result, &error); + + // printf ("result = % .18f\n", result); + // printf ("estimated error = % .18f\n", error); + // printf ("intervals = %zu\n", w->size); + + gsl_integration_workspace_free (w); + + return fabs(result); +} + +/* Calculate the total attenuation for each type of ray */ +double IceRayTracing::GetTotalAttenuationDirect (double A0, double frequency, double z0, double z1, double Lvalue) { + z0=fabs(z0); + z1=fabs(z1); + return IceRayTracing::IntegrateOverLAttn(A0,frequency,z0,z1,Lvalue); +} + +double IceRayTracing::GetTotalAttenuationReflected (double A0, double frequency, double z0, double z1, double Lvalue) { + z0=fabs(z0); + z1=fabs(z1); + return IceRayTracing::IntegrateOverLAttn(A0,frequency,z0,0.000001,Lvalue) + IceRayTracing::IntegrateOverLAttn(A0,frequency,z1,0.000001,Lvalue); +} + +double IceRayTracing::GetTotalAttenuationRefracted (double A0, double frequency, double z0, double z1, double zmax, double Lvalue) { + z0=fabs(z0); + z1=fabs(z1); + return IceRayTracing::IntegrateOverLAttn(A0,frequency,z0,zmax,Lvalue) + IceRayTracing::IntegrateOverLAttn(A0,frequency,z1,zmax,Lvalue); +} + +/* Use GSL minimiser which relies on calculating function deriavtives. This function uses GSL's Newton's algorithm to find root for a given function. */ +double IceRayTracing::FindFunctionRootFDF(gsl_function_fdf FDF,double x_lo, double x_hi){ + int status; + int iter = 0, max_iter = 100; + const gsl_root_fdfsolver_type *T; + gsl_root_fdfsolver *s; + double x0 ,x = (x_lo+x_hi)/2; + + T = gsl_root_fdfsolver_newton; + s = gsl_root_fdfsolver_alloc (T); + gsl_root_fdfsolver_set (s, &FDF, x); + + // printf ("using %s method\n", + // gsl_root_fdfsolver_name (s)); + + // printf ("%-5s %10s %10s %10s\n", + // "iter", "root", "err", "err(est)"); + do + { + iter++; + status = gsl_root_fdfsolver_iterate (s); + x0 = x; + x = gsl_root_fdfsolver_root (s); + status = gsl_root_test_delta (x, x0, 0, 1e-6); + + // if (status == GSL_SUCCESS) + // printf ("Converged:\n"); + + // printf ("%5d %10.7f %10.7f\n", + // iter, x, x - x0); + + } + while (status == GSL_CONTINUE && iter < max_iter); + + gsl_root_fdfsolver_free (s); + + return x; +} + +/* Use GSL minimiser which uses GSL's false position algorithm to find root for a given function. */ +double IceRayTracing::FindFunctionRoot(gsl_function F,double x_lo, double x_hi) +{ + int status; + int iter = 0, max_iter = 100; + const gsl_root_fsolver_type *T; + gsl_root_fsolver *s; + double r = 0; + + T = gsl_root_fsolver_falsepos; + s = gsl_root_fsolver_alloc (T); + gsl_set_error_handler_off(); + status = gsl_root_fsolver_set (s, &F, x_lo, x_hi); + + do + { + iter++; + status = gsl_root_fsolver_iterate (s); + r = gsl_root_fsolver_root (s); + x_lo = gsl_root_fsolver_x_lower (s); + x_hi = gsl_root_fsolver_x_upper (s); + + if(x_lo>x_hi){ + swap(x_lo,x_hi); + } + double checkval=(*((F).function))(r,(F).params); + status = gsl_root_test_residual (checkval,1e-6); + //status = gsl_root_test_interval (x_lo, x_hi,0, 0.000001); + + if (status == GSL_SUCCESS){ + // printf ("Converged:"); + // printf ("%5d [%.7f, %.7f] %.7f %.7f\n",iter, x_lo, x_hi,r,x_hi - x_lo); + } + } + while (status == GSL_CONTINUE && iter < max_iter); + + //printf ("status = %s\n", gsl_strerror (status)); + gsl_root_fsolver_free (s); + + return r; +} + +/* Use GSL minimiser which uses GSL's false position algorithm to find root for a given function. */ +double IceRayTracing::FindFunctionRootZmax(gsl_function F,double x_lo, double x_hi) +{ + int status; + int iter = 0, max_iter = 100; + const gsl_root_fsolver_type *T; + gsl_root_fsolver *s; + double r = 0; + + T = gsl_root_fsolver_falsepos; + s = gsl_root_fsolver_alloc (T); + gsl_set_error_handler_off(); + status = gsl_root_fsolver_set (s, &F, x_lo, x_hi); + + do + { + iter++; + status = gsl_root_fsolver_iterate (s); + r = gsl_root_fsolver_root (s); + x_lo = gsl_root_fsolver_x_lower (s); + x_hi = gsl_root_fsolver_x_upper (s); + status = gsl_root_test_interval (x_lo, x_hi,1e-6, 1e-6); + + if (status == GSL_SUCCESS){ + // printf ("Converged:"); + // printf ("%5d [%.7f, %.7f] %.7f %.7f\n",iter, x_lo, x_hi,r,x_hi - x_lo); + } + } + while (status == GSL_CONTINUE && iter < max_iter); + + gsl_root_fsolver_free (s); + + return r; +} + +/* Define the function that will be minimized to get the value of the depth of the turning point for a given refracted ray. This function basically requires the value of the L parameter to find the depth of the turning point. This comes from the part of the fDnfR function where sqrt( n(z) - L ). This imposes the constraint then n(z)=L at the turning point of the ray from which we can find zmax. */ +double IceRayTracing::GetMinnz(double x,void *params){ + struct IceRayTracing::Minnz_params *p= (struct IceRayTracing::Minnz_params *) params; + double A = p->a; + double L = p->l; + return A+IceRayTracing::GetB(x)*exp(-IceRayTracing::GetC(x)*x)-L; +} + +/* Get the value of the depth of the turning point for the refracted ray */ +double IceRayTracing::GetZmax(double A, double L){ + gsl_function F1; + struct IceRayTracing::Minnz_params params1= {A,L}; + F1.function = &GetMinnz; + F1.params = ¶ms1; + double zmax=IceRayTracing::FindFunctionRootZmax(F1,0.0,5000); + return zmax; +} + +/* Analytical solution describing ray paths in ice as function of depth */ +double IceRayTracing::fDnfR(double x,void *params){ + + struct IceRayTracing::fDnfR_params *p= (struct IceRayTracing::fDnfR_params *) params; + double A = p->a; + double B = p->b; + double C = p->c; + double L = p->l; + + return (L/C)*(1.0/sqrt(A*A-L*L))*(C*x-log(A*IceRayTracing::Getnz(x)-L*L+sqrt(A*A-L*L)*sqrt(pow(IceRayTracing::Getnz(x),2)-L*L)));; +} + +/* Analytical solution describing the ray path in ice as a function of the L parameter */ +double IceRayTracing::fDnfR_L(double x,void *params){ + + struct IceRayTracing::fDnfR_L_params *p= (struct IceRayTracing::fDnfR_L_params *) params; + double A = p->a; + double B = p->b; + double C = p->c; + double Z = p->z; + + double result=(x/C)*(1.0/sqrt(A*A-x*x))*(C*Z-log(A*IceRayTracing::Getnz(Z)-x*x+sqrt(A*A-x*x)*sqrt(pow(IceRayTracing::Getnz(Z),2)-x*x))); + + return result; +} + +/* The function used to calculate ray propogation time in ice */ +double IceRayTracing::ftimeD(double x,void *params){ + + struct IceRayTracing::ftimeD_params *p= (struct IceRayTracing::ftimeD_params *) params; + double A = p->a; + double B = p->b; + double C = p->c; + double Speedc = p->speedc; + double L = p->l; + + return (1.0/(Speedc*C*sqrt(pow(IceRayTracing::Getnz(x),2)-L*L)))*(pow(IceRayTracing::Getnz(x),2)-L*L+(C*x-log(A*IceRayTracing::Getnz(x)-L*L+sqrt(A*A-L*L)*sqrt(pow(IceRayTracing::Getnz(x),2)-L*L)))*(A*A*sqrt(pow(IceRayTracing::Getnz(x),2)-L*L))/sqrt(A*A-L*L) +A*sqrt(pow(IceRayTracing::Getnz(x),2)-L*L)*log(IceRayTracing::Getnz(x)+sqrt(pow(IceRayTracing::Getnz(x),2)-L*L)) ); +} + +/* The function is used to calculate ray geometric path in ice */ +double IceRayTracing::fpathD(double x,void *params){ + + struct IceRayTracing::ftimeD_params *p= (struct IceRayTracing::ftimeD_params *) params; + double A = p->a; + double B = p->b; + double C = p->c; + double Speedc = p->speedc; + double L = p->l; + + //integral sec(sin^(-1)(L/(A + B e^(C x)))) dx = (log((A + B e^(C x)) (sqrt((A^2 + 2 A B e^(C x) + B^2 e^(2 C x) - L^2)/(A + B e^(C x))^2) + 1)) - (A log(A sqrt(A^2 - L^2) sqrt((A^2 + 2 A B e^(C x) + B^2 e^(2 C x) - L^2)/(A + B e^(C x))^2) + B sqrt(A^2 - L^2) e^(C x) sqrt((A^2 + 2 A B e^(C x) + B^2 e^(2 C x) - L^2)/(A + B e^(C x))^2) + A^2 + A B e^(C x) - L^2))/sqrt(A^2 - L^2) + (A C x)/sqrt(A^2 - L^2))/C; + + return (log((A + B*exp(C*x))*(sqrt((A*A + 2*A*B*exp(C*x) + B*B*exp(2*C*x) - L*L)/((A + B*exp(C*x))*(A + B*exp(C*x))) ) + 1)) - (A*log(A*sqrt(A*A - L*L)*sqrt((A*A + 2*A*B*exp(C*x) + B*B* exp(2*C*x) - L*L)/(( A + B*exp(C*x))*(A + B*exp(C*x)))) + B*sqrt(A*A - L*L)*exp(C*x)*sqrt((A*A + 2*A*B*exp(C*x) + B*B* exp(2*C*x) - L*L)/((A + B*exp(C*x))*(A + B*exp(C*x)))) + A*A + A*B*exp(C*x) - L*L))/sqrt(A*A - L*L) + (A*C*x)/sqrt(A*A - L*L))/C ; + +} + +/* The set of functions starting with the name "fDa" are used in the minimisation procedure to find the launch angle (or the L parameter) for the direct ray */ +double IceRayTracing::fDa(double x,void *params){ + struct IceRayTracing::fDanfRa_params *p= (struct IceRayTracing::fDanfRa_params *) params; + double A = p->a; + double z0 = p->z0; + double x1 = p->x1; + double z1 = p->z1; + + struct IceRayTracing::fDnfR_L_params params1a = {A, IceRayTracing::GetB(z1), IceRayTracing::GetC(z1), z1}; + struct IceRayTracing::fDnfR_L_params params1b = {A, IceRayTracing::GetB(z0), IceRayTracing::GetC(z0), z0}; + struct IceRayTracing::fDnfR_L_params params1c = {A, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), IceRayTracing::GetC(IceRayTracing::TransitionBoundary), -IceRayTracing::TransitionBoundary}; + struct IceRayTracing::fDnfR_L_params params1d = {A, IceRayTracing::GetB(-(IceRayTracing::TransitionBoundary+0.000001)), IceRayTracing::GetC(-(IceRayTracing::TransitionBoundary+0.000001)), -(IceRayTracing::TransitionBoundary+0.000001)}; + + double distancez0z1=0; + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)>IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)IceRayTracing::TransitionBoundary && fabs(z1)==IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1c) + IceRayTracing::fDnfR_L(x,¶ms1d) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + }else{ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + // if(std::isnan(distancez0z1)){ + // distancez0z1=1e9; + // } + double output=distancez0z1-x1; + + return output; +} + +double IceRayTracing::fDa_df(double x,void *params){ + gsl_function F; + F.function = &IceRayTracing::fDa; + F.params = params; + + double result,abserr; + gsl_deriv_central (&F, x, 1e-8, &result, &abserr); + + return result; +} + +void IceRayTracing::fDa_fdf (double x, void *params,double *y, double *dy){ + *y = IceRayTracing::fDa(x,params); + *dy = IceRayTracing::fDa_df(x,params); +} + +/* The set of functions starting with the name "fRa" are used in the minimisation procedure to find the launch angle (or the L parameter) for the reflected ray */ +double IceRayTracing::fRa(double x,void *params){ + struct IceRayTracing::fDanfRa_params *p= (struct IceRayTracing::fDanfRa_params *) params; + double A = p->a; + double z0 = p->z0; + double x1 = p->x1; + double z1 = p->z1; + + struct IceRayTracing::fDnfR_L_params params1a = {A, IceRayTracing::GetB(z1), -IceRayTracing::GetC(z1), -z1}; + struct IceRayTracing::fDnfR_L_params params1b = {A, IceRayTracing::GetB(z0), -IceRayTracing::GetC(z0), -z0}; + struct IceRayTracing::fDnfR_L_params params1c = {A, IceRayTracing::GetB(1e-7), -IceRayTracing::GetC(1e-7), 1e-7}; + struct IceRayTracing::fDnfR_L_params params1d = {A, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary), IceRayTracing::TransitionBoundary}; + struct IceRayTracing::fDnfR_L_params params1f = {A, IceRayTracing::GetB(IceRayTracing::TransitionBoundary+0.000001), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary+0.000001), IceRayTracing::TransitionBoundary+0.000001}; + + double distancez0z1=0; + double distancez0surface=0; + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)>IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1d) + IceRayTracing::fDnfR_L(x,¶ms1f) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)IceRayTracing::TransitionBoundary && fabs(z1)==IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1d) + IceRayTracing::fDnfR_L(x,¶ms1f) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1d) + IceRayTracing::fDnfR_L(x,¶ms1f) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + }else{ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + // if(std::isnan(distancez0z1)){ + // distancez0z1=1e9; + // } + // if(std::isnan(distancez0surface)){ + // distancez0surface=1e9; + // } + double output= distancez0z1 - 2*(distancez0surface) - x1; + + return output; +} + +double IceRayTracing::fRa_df(double x,void *params){ + gsl_function F; + F.function = &IceRayTracing::fRa; + F.params = params; + + double result,abserr; + gsl_deriv_central (&F, x, 1e-8, &result, &abserr); + + return result; +} + +void IceRayTracing::fRa_fdf (double x, void *params,double *y, double *dy){ + *y = IceRayTracing::fRa(x,params); + *dy = IceRayTracing::fRa_df(x,params); +} + +/* This function is minimised to find the launch angle (or the L parameter) for the refracted ray */ +double IceRayTracing::fRaa(double x,void *params){ + struct IceRayTracing::fDanfRa_params *p= (struct IceRayTracing::fDanfRa_params *) params; + double A = p->a; + double z0 = p->z0; + double x1 = p->x1; + double z1 = p->z1; + + double zmax= IceRayTracing::GetZmax(A,x)+1e-7; + double output=0; + if(zmax>0){ + struct IceRayTracing::fDnfR_L_params params1a = {A, IceRayTracing::GetB(z1), -IceRayTracing::GetC(z1), -z1}; + struct IceRayTracing::fDnfR_L_params params1b = {A, IceRayTracing::GetB(z0), -IceRayTracing::GetC(z0), -z0}; + struct IceRayTracing::fDnfR_L_params params1c = {A, IceRayTracing::GetB(zmax), -IceRayTracing::GetC(zmax), zmax}; + struct IceRayTracing::fDnfR_L_params params1d = {A, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary), IceRayTracing::TransitionBoundary}; + struct IceRayTracing::fDnfR_L_params params1f = {A, IceRayTracing::GetB(IceRayTracing::TransitionBoundary+1e-7), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary+1e-7), IceRayTracing::TransitionBoundary+1e-7}; + + double distancez0z1=0; + double distancez0surface=0; + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)>IceRayTracing::TransitionBoundary){ + if(zmax<=IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1d) + IceRayTracing::fDnfR_L(x,¶ms1f) - IceRayTracing::fDnfR_L(x,¶ms1b); + }else{ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)IceRayTracing::TransitionBoundary && fabs(z1)==IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1d) + IceRayTracing::fDnfR_L(x,¶ms1f) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1d) + IceRayTracing::fDnfR_L(x,¶ms1f) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + }else{ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + distancez0surface=IceRayTracing::fDnfR_L(x,¶ms1c) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + + if(std::isnan(distancez0z1)){ + distancez0z1=1e9; + } + if(std::isnan(distancez0surface)){ + distancez0surface=1e9; + } + output= distancez0z1 - 2*(distancez0surface) - x1; + }else{ + output=1e9; + } + return output; +} + +double IceRayTracing::fRaa_df(double x,void *params){ + gsl_function F; + F.function = &IceRayTracing::fRaa; + F.params = params; + + double result,abserr; + gsl_deriv_central (&F, x, 1e-8, &result, &abserr); + + return result; +} + +void IceRayTracing::fRaa_fdf (double x, void *params,double *y, double *dy){ + *y = IceRayTracing::fRaa(x,params); + *dy = IceRayTracing::fRaa_df(x,params); +} + +/* This functions works for the Direct ray and gives you back the launch angle, receive angle and propagation time of the ray together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. */ +double* IceRayTracing::GetDirectRayPar(double z0, double x1, double z1){ + + double *output=new double[6]; + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* First we setup the fDa function that will be minimised to get the launch angle (or the L parameter) for the direct ray. */ + gsl_function F1; + struct IceRayTracing::fDanfRa_params params1= {IceRayTracing::A_ice, z0, x1, z1}; + F1.function = &IceRayTracing::fDa; + F1.params = ¶ms1; + + // gsl_function_fdf F1; + // struct IceRayTracing::fDanfRa_params params1= {IceRayTracing::A_ice, z0, x1, z1}; + // F1.f = &fDa; + // F1.df = &fDa_df; + // F1.fdf = &fDa_fdf; + // F1.params = ¶ms1; + + /* In my raytracing solution given in the function IceRayTracing::fDnfR the launch angle (or the L parameter) has limit placed on it by this part in the solution sqrt( n(z)^2 - L^2) . This sqrt cannot be negative for both z0 and z1 and this sets the upper limit in our minimisation to get the launch angle (or the L parameter). Here I am basically setting the upper limit as GSL requires that my function is well behaved on the upper and lower bounds I give it for minimisation. */ + double UpLimnz[]={IceRayTracing::Getnz(z1),IceRayTracing::Getnz(z0)}; + double* UpperLimitL=min_element(UpLimnz,UpLimnz+2); + + /* Do the minimisation and get the value of the L parameter and the launch angle and then verify to see that the value of L that we got was actually a root of fDa function. */ + double lvalueD=IceRayTracing::FindFunctionRoot(F1,1e-7,UpperLimitL[0]); + double LangD=asin(lvalueD/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + double checkzeroD=IceRayTracing::fDa(lvalueD,¶ms1); + + /* Get the propagation time for the direct ray using the ftimeD function after we have gotten the value of the L parameter. */ + struct IceRayTracing::ftimeD_params params2a = {IceRayTracing::A_ice, IceRayTracing::GetB(z0), -IceRayTracing::GetC(z0), IceRayTracing::c_light_ms,lvalueD}; + struct IceRayTracing::ftimeD_params params2b = {IceRayTracing::A_ice, IceRayTracing::GetB(z1), -IceRayTracing::GetC(z1), IceRayTracing::c_light_ms,lvalueD}; + struct IceRayTracing::ftimeD_params params2c = {IceRayTracing::A_ice, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary), IceRayTracing::c_light_ms, lvalueD}; + struct IceRayTracing::ftimeD_params params2d = {IceRayTracing::A_ice, IceRayTracing::GetB(IceRayTracing::TransitionBoundary+1e-7), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary+1e-7), IceRayTracing::c_light_ms, lvalueD}; + + /* we do the subtraction because we are measuring the time taken between the Tx and Rx positions */ + double timeD=0; + double pathD=0; + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)>IceRayTracing::TransitionBoundary){ + timeD=IceRayTracing::ftimeD(-z0,¶ms2a) - IceRayTracing::ftimeD(-z1,¶ms2b); + pathD=IceRayTracing::fpathD(-z0,¶ms2a) - IceRayTracing::fpathD(-z1,¶ms2b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)IceRayTracing::TransitionBoundary && fabs(z1)==IceRayTracing::TransitionBoundary){ + timeD=IceRayTracing::ftimeD(-z0,¶ms2a) - IceRayTracing::ftimeD(IceRayTracing::TransitionBoundary+1e-7,¶ms2d) + IceRayTracing::ftimeD(IceRayTracing::TransitionBoundary,¶ms2c) - IceRayTracing::ftimeD(-z1,¶ms2b); + pathD=IceRayTracing::fpathD(-z0,¶ms2a) - IceRayTracing::fpathD(IceRayTracing::TransitionBoundary+1e-7,¶ms2d) + IceRayTracing::fpathD(IceRayTracing::TransitionBoundary,¶ms2c) - IceRayTracing::fpathD(-z1,¶ms2b); + } + }else{ + timeD=IceRayTracing::ftimeD(-z0,¶ms2a) - IceRayTracing::ftimeD(-z1,¶ms2b); + pathD=IceRayTracing::fpathD(-z0,¶ms2a) - IceRayTracing::fpathD(-z1,¶ms2b); + } + + /* Setup the function that will be used to calculate the angle of reception for all the rays */ + gsl_function F5; + struct IceRayTracing::fDnfR_params params5a = {IceRayTracing::A_ice, IceRayTracing::GetB(z1), -IceRayTracing::GetC(z1), lvalueD}; + double result, abserr; + F5.function = &IceRayTracing::fDnfR; + + /* Calculate the recieve angle for direc rays by calculating the derivative of the function at the Rx position */ + F5.params = ¶ms5a; + gsl_deriv_central (&F5, -z1, 1e-8, &result, &abserr); + double RangD=atan(result)*(180.0/IceRayTracing::pi); + + /* When the Tx and Rx are at the same depth my function struggles to find a ray between them when they are very close to each other. In that case the ray is pretty much like a straight line. */ + if(z1==z0 && std::isnan(RangD)==true){ + RangD=180-LangD; + } + + /* This sometimes happens that when the Rx is very close to the peak point (or the turning point) of the ray then its hard to calculate the derivative around that area since the solution blows up around that area. therefore this is a good approximation. */ + if(z1!=z0 && std::isnan(RangD)==true){ + RangD=90; + } + + dsw=0; + /* If the Tx and Rx depth were switched then put them back to their original position */ + if(Flip==true){ + dsw=z0; + z0=z1; + z1=dsw; + } + + output[0]=RangD; + output[1]=LangD; + output[2]=timeD; + output[3]=lvalueD; + output[4]=checkzeroD; + output[5]=pathD; + + /* If the flip case is true where we flipped Rx and Tx depths to trace rays then make sure everything is switched back before we give the output to the user. */ + if(Flip==true){ + output[0]=180-LangD; + output[1]=180-RangD; + } + + return output; +} + +/* This functions works for the Reflected ray and gives you back the launch angle, receive angle and propagation times (of the whole ray and the two direct rays that make it up) together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. */ +double *IceRayTracing::GetReflectedRayPar(double z0, double x1 ,double z1){ + + double *output=new double[11]; + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* First we setup the fRa function that will be minimised to get the launch angle (or the L parameter) for the reflected ray. */ + gsl_function F3; + struct IceRayTracing::fDanfRa_params params3= {IceRayTracing::A_ice, z0, x1, z1}; + F3.function = &IceRayTracing::fRa; + F3.params = ¶ms3; + + // gsl_function_fdf F3; + // struct IceRayTracing::fDanfRa_params params3= {IceRayTracing::A_ice, z0, x1, z1}; + // F3.f = &IceRayTracing::fRa; + // F3.df = &IceRayTracing::fRa_df; + // F3.fdf = &IceRayTracing::fRa_fdf; + // F3.params = ¶ms3; + + /* In my raytracing solution given in the function IceRayTracing::fDnfR the launch angle (or the L parameter) has limit placed on it by this part in the solution sqrt( n(z)^2 - L^2) . This sqrt cannot be negative for both z0, z1 and also 1e-7 m and this sets the upper limit in our minimisation to get the launch angle (or the L parameter). Here I am basically setting the upper limit as GSL requires that my function is well behaved on the upper and lower bounds I give it for minimisation. */ + double UpLimnz[]={IceRayTracing::Getnz(z1),IceRayTracing::Getnz(z0),IceRayTracing::Getnz(1e-7)}; + double *UpperLimitL=min_element(UpLimnz,UpLimnz+3); + + /* Do the minimisation and get the value of the L parameter and the launch angle and then verify to see that the value of L that we got was actually a root of fRa function. */ + double lvalueR=IceRayTracing::FindFunctionRoot(F3,1e-7,UpperLimitL[0]); + double LangR=asin(lvalueR/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + double checkzeroR=IceRayTracing::fRa(lvalueR,¶ms3); + + /* Get the propagation time for the reflected ray using the ftimeD function after we have gotten the value of the L parameter. */ + struct IceRayTracing::ftimeD_params params3a = {IceRayTracing::A_ice, IceRayTracing::GetB(z0), IceRayTracing::GetC(z0), IceRayTracing::c_light_ms,lvalueR}; + struct IceRayTracing::ftimeD_params params3b = {IceRayTracing::A_ice, IceRayTracing::GetB(z1), IceRayTracing::GetC(z1), IceRayTracing::c_light_ms,lvalueR}; + struct IceRayTracing::ftimeD_params params3c = {IceRayTracing::A_ice, IceRayTracing::GetB(1e-7), IceRayTracing::GetC(1e-7), IceRayTracing::c_light_ms,lvalueR}; + struct IceRayTracing::ftimeD_params params3d = {IceRayTracing::A_ice, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), IceRayTracing::GetC(IceRayTracing::TransitionBoundary), IceRayTracing::c_light_ms, lvalueR}; + struct IceRayTracing::ftimeD_params params3f = {IceRayTracing::A_ice, IceRayTracing::GetB(IceRayTracing::TransitionBoundary+1e-7), IceRayTracing::GetC(IceRayTracing::TransitionBoundary+1e-7), IceRayTracing::c_light_ms, lvalueR}; + /* We do the subtraction because we are measuring the time taken between the Tx and Rx positions. In the reflected case we basically have two direct rays 1) from Tx to surface 2) from surface to Rx. Also get the time for the two individual direct rays separately */ + double timeR1=0; + double timeR2=0; + double pathR1=0; + double pathR2=0; + + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)>IceRayTracing::TransitionBoundary){ + timeR1=IceRayTracing::ftimeD(-1e-7,¶ms3c) - IceRayTracing::ftimeD(-IceRayTracing::TransitionBoundary,¶ms3d) + IceRayTracing::ftimeD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms3f) - IceRayTracing::ftimeD(z0,¶ms3a); + timeR2=IceRayTracing::ftimeD(-1e-7,¶ms3c) - IceRayTracing::ftimeD(-IceRayTracing::TransitionBoundary,¶ms3d) + IceRayTracing::ftimeD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms3f) - IceRayTracing::ftimeD(z1,¶ms3b); + + pathR1=IceRayTracing::fpathD(-1e-7,¶ms3c) - IceRayTracing::fpathD(-IceRayTracing::TransitionBoundary,¶ms3d) + IceRayTracing::fpathD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms3f) - IceRayTracing::fpathD(z0,¶ms3a); + pathR2=IceRayTracing::fpathD(-1e-7,¶ms3c) - IceRayTracing::fpathD(-IceRayTracing::TransitionBoundary,¶ms3d) + IceRayTracing::fpathD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms3f) - IceRayTracing::fpathD(z1,¶ms3b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)IceRayTracing::TransitionBoundary && fabs(z1)==IceRayTracing::TransitionBoundary){ + timeR1=IceRayTracing::ftimeD(-1e-7,¶ms3c) - IceRayTracing::ftimeD(-IceRayTracing::TransitionBoundary,¶ms3d) + IceRayTracing::ftimeD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms3f) - IceRayTracing::ftimeD(z0,¶ms3a); + timeR2=IceRayTracing::ftimeD(-1e-7,¶ms3c) - IceRayTracing::ftimeD(z1,¶ms3b); + + pathR1=IceRayTracing::fpathD(-1e-7,¶ms3c) - IceRayTracing::fpathD(-IceRayTracing::TransitionBoundary,¶ms3d) + IceRayTracing::fpathD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms3f) - IceRayTracing::fpathD(z0,¶ms3a); + pathR2=IceRayTracing::fpathD(-1e-7,¶ms3c) - IceRayTracing::fpathD(z1,¶ms3b); + } + }else{ + timeR1=IceRayTracing::ftimeD(-1e-7,¶ms3c) - IceRayTracing::ftimeD(z0,¶ms3a); + timeR2=IceRayTracing::ftimeD(-1e-7,¶ms3c) - IceRayTracing::ftimeD(z1,¶ms3b); + + pathR1=IceRayTracing::fpathD(-1e-7,¶ms3c) - IceRayTracing::fpathD(z0,¶ms3a); + pathR2=IceRayTracing::fpathD(-1e-7,¶ms3c) - IceRayTracing::fpathD(z1,¶ms3b); + } + double timeR=timeR1 + timeR2; + double pathR=pathR1 + pathR2; + + /* flip the times back if the original positions were flipped */ + if(Flip==true){ + double dumR=timeR2; + timeR2=timeR1; + timeR1=dumR; + + double dumRb=pathR2; + pathR2=pathR1; + pathR1=dumRb; + } + timeR1=timeR1; + timeR2=timeR2; + + pathR1=pathR1; + pathR2=pathR2; + + /* Setup the function that will be used to calculate the angle of reception for all the rays */ + gsl_function F5; + struct IceRayTracing::fDnfR_params params5b = {IceRayTracing::A_ice, IceRayTracing::GetB(z1), IceRayTracing::GetC(z1), lvalueR}; + double result, abserr; + F5.function = &IceRayTracing::fDnfR; + + /* Calculate the recieve angle for reflected ray by calculating the derivative of the function at the Rx position */ + F5.params = ¶ms5b; + gsl_deriv_central (&F5, z1, 1e-8, &result, &abserr); + double RangR=180-atan(result)*(180.0/IceRayTracing::pi); + + /* When the Tx and Rx are at the same depth my function struggles to find a ray between them when they are very close to each other. In that case the ray is pretty much like a straight line. */ + if(z1==z0 && std::isnan(RangR)==true){ + RangR=180-LangR; + } + + /* This sometimes happens that when the Rx is very close to the peak point (or the turning point) of the ray then its hard to calculate the derivative around that area since the solution blows up around that area. therefore this is a good approximation. */ + if(z1!=z0 && std::isnan(RangR)==true){ + RangR=90; + } + + if(std::isnan(checkzeroR)){ + checkzeroR=-1000; + } + + /* Calculate the angle of incidence of the reflected ray at the surface ice. This will be used to calculate the Fresnel Coefficients. The angle is calculated by calculating the derivative of the ray path fucnction at the surface*/ + struct IceRayTracing::fDnfR_params paramsIAngB = {IceRayTracing::A_ice, IceRayTracing::GetB(1e-7), IceRayTracing::GetC(1e-7), lvalueR}; + F5.function = &IceRayTracing::fDnfR; + F5.params = ¶msIAngB; + gsl_deriv_central (&F5, -1e-7, 1e-8, &result, &abserr); + double IncidenceAngleInIce=atan(result)*(180.0/IceRayTracing::pi); + + dsw=0; + /* If the Tx and Rx depth were switched then put them back to their original position */ + if(Flip==true){ + dsw=z0; + z0=z1; + z1=dsw; + } + + output[0]=RangR; + output[1]=LangR; + output[2]=timeR; + output[3]=lvalueR; + output[4]=checkzeroR; + output[5]=timeR1; + output[6]=timeR2; + output[7]=IncidenceAngleInIce; + output[8]=pathR; + output[9]=pathR1; + output[10]=pathR2; + + /* If the flip case is true where we flipped Rx and Tx depths to trace rays then make sure everything is switched back before we give the output to the user. */ + if(Flip==true){ + output[0]=180-LangR; + output[1]=180-RangR; + } + + return output; +} + +/* This functions works for the Refracted ray and gives you back the launch angle, receive angle and propagation times (of the whole ray and the two direct rays that make it up) together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. It requires the launch angle of the reflected ray as an input. */ +double *IceRayTracing::GetRefractedRayPar(double z0, double x1 ,double z1, double LangR, double RangR, double checkzeroD, double checkzeroR){ + + double *output=new double[22]; + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + dsw=0; + if(Flip==true){ + dsw=180-LangR; + LangR=180-RangR; + RangR=dsw; + } + + /* Set up all the variables that will be used to get the parameters for the refracted ray */ + //double lvalueR=sin(LangR*(IceRayTracing::pi/180))*IceRayTracing::Getnz(z0); + double lvalueRa[2]={0,0}; + double LangRa[2]={0,0}; + double checkzeroRa[2]={-1000,-1000}; + + double timeRa[2]={0,0}; + double timeRa1[2]={0,0}; + double timeRa2[2]={0,0}; + + double pathRa[2]={0,0}; + double pathRa1[2]={0,0}; + double pathRa2[2]={0,0}; + + double raytime[2]={0,0}; + double RangRa[2]={0,0}; + double zmax[2]={10,10}; + + /* In my raytracing solution given in the function IceRayTracing::fDnfR the launch angle (or the L parameter) has limit placed on it by this part in the solution sqrt( n(z)^2 - L^2) . This sqrt cannot be negative for both z0, z1 and also 1e-7 m and this sets the upper limit in our minimisation to get the launch angle (or the L parameter). Here I am basically setting the upper limit as GSL requires that my function is well behaved on the upper and lower bounds I give it for minimisation. */ + double UpLimnz[]={IceRayTracing::Getnz(z0),IceRayTracing::Getnz(z1)}; + double *UpperLimitL=min_element(UpLimnz,UpLimnz+2); + + /* First we setup the fRa function that will be minimised to get the launch angle (or the L parameter) for the refracted ray. */ + gsl_function F4; + struct IceRayTracing::fDanfRa_params params4= {IceRayTracing::A_ice, z0, x1, z1}; + F4.function = &IceRayTracing::fRaa; + F4.params = ¶ms4; + + gsl_function_fdf F4b; + F4b.f = &IceRayTracing::fRaa; + F4b.df = &IceRayTracing::fRaa_df; + F4b.fdf = &IceRayTracing::fRaa_fdf; + F4b.params = ¶ms4; + + /* Do the minimisation and get the value of the L parameter and the launch angle and then verify to see that the value of L that we got was actually a root of fRaa function. The thing to note here is the lower limit of the minimisation function is set to the L value corresponding to the reflected ray launch angle. Since we know the refracted ray always has bigger launch angle the reflected ray this reduces our range and makes the function more efficient at finding the refracted ray launch angle. */ + double LowerLimit=0; + LowerLimit=IceRayTracing::Getnz(z0)*sin((64.0*(IceRayTracing::pi/180.0))); + if(LowerLimit>UpperLimitL[0]){ + LowerLimit=IceRayTracing::Getnz(z0)*sin((LangR*(IceRayTracing::pi/180.0))); + } + + lvalueRa[0]=IceRayTracing::FindFunctionRoot(F4,LowerLimit,UpperLimitL[0]); + LangRa[0]=asin(lvalueRa[0]/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + checkzeroRa[0]=(IceRayTracing::fRaa(lvalueRa[0],¶ms4)); + zmax[0]=IceRayTracing::GetZmax(IceRayTracing::A_ice,lvalueRa[0])+1e-7; + + if(fabs(checkzeroRa[0])>0.5){ + lvalueRa[0]=IceRayTracing::FindFunctionRootFDF(F4b,LowerLimit,UpperLimitL[0]); + LangRa[0]=asin(lvalueRa[0]/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + checkzeroRa[0]=IceRayTracing::fRaa(lvalueRa[0],¶ms4); + zmax[0]=IceRayTracing::GetZmax(IceRayTracing::A_ice,lvalueRa[0])+1e-7; + } + + if(lvalueRa[0]<0){ + checkzeroRa[0]=-1000; + } + + if(fabs(checkzeroRa[0])<0.5 && fabs(checkzeroD)>0.5 && fabs(checkzeroR)>0.5){ + lvalueRa[1]=IceRayTracing::FindFunctionRoot(F4,lvalueRa[0]-0.23,lvalueRa[0]-0.023); + LangRa[1]=asin(lvalueRa[1]/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + checkzeroRa[1]=IceRayTracing::fRaa(lvalueRa[1],¶ms4); + zmax[1]=IceRayTracing::GetZmax(IceRayTracing::A_ice,lvalueRa[1])+1e-7; + + if(fabs(checkzeroRa[1])>0.5 || std::isnan(checkzeroRa[1])==true || fabs(lvalueRa[1]-lvalueRa[0])<1e-4 ){ + lvalueRa[1]=IceRayTracing::FindFunctionRoot(F4,lvalueRa[0]-0.15,lvalueRa[0]-0.023); + LangRa[1]=asin(lvalueRa[1]/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + checkzeroRa[1]=IceRayTracing::fRaa(lvalueRa[1],¶ms4); + zmax[1]=IceRayTracing::GetZmax(IceRayTracing::A_ice,lvalueRa[1])+1e-7; + } + + if(fabs(checkzeroRa[1])>0.5 || std::isnan(checkzeroRa[1])==true || fabs(lvalueRa[1]-lvalueRa[0])<1e-4 ){ + if(lvalueRa[0]+0.0050.5 || std::isnan(checkzeroRa[1])==true || fabs(lvalueRa[1]-lvalueRa[0])<1e-4 ){ + double lvalueRatmp=IceRayTracing::FindFunctionRootFDF(F4b,lvalueRa[0]-0.23,lvalueRa[0]-0.023); + if(fabs(lvalueRatmp)0.5 || std::isnan(checkzeroRa[1])==true || fabs(lvalueRa[1]-lvalueRa[0])<1e-4 ){ + double lvalueRatmp=IceRayTracing::FindFunctionRootFDF(F4b,lvalueRa[0]-0.1,lvalueRa[0]-0.023); + if(fabs(lvalueRatmp)IceRayTracing::TransitionBoundary && fabs(z1)>IceRayTracing::TransitionBoundary){ + if(zmax[i]<=IceRayTracing::TransitionBoundary){ + timeRa1[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(-IceRayTracing::TransitionBoundary,¶ms4d) + IceRayTracing::ftimeD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms4f) - IceRayTracing::ftimeD(z0,¶ms4a); + timeRa2[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(-IceRayTracing::TransitionBoundary,¶ms4d) + IceRayTracing::ftimeD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms4f) - IceRayTracing::ftimeD(z1,¶ms4b); + + pathRa1[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(-IceRayTracing::TransitionBoundary,¶ms4d) + IceRayTracing::fpathD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms4f) - IceRayTracing::fpathD(z0,¶ms4a); + pathRa2[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(-IceRayTracing::TransitionBoundary,¶ms4d) + IceRayTracing::fpathD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms4f) - IceRayTracing::fpathD(z1,¶ms4b); + + }else{ + timeRa1[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(z0,¶ms4a); + timeRa2[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(z1,¶ms4b); + + pathRa1[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(z0,¶ms4a); + pathRa2[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(z1,¶ms4b); + } + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(z1)IceRayTracing::TransitionBoundary && fabs(z1)==IceRayTracing::TransitionBoundary){ + timeRa1[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(-IceRayTracing::TransitionBoundary,¶ms4d) + IceRayTracing::ftimeD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms4f) - IceRayTracing::ftimeD(z0,¶ms4a); + timeRa2[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(z1,¶ms4b); + + pathRa1[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(-IceRayTracing::TransitionBoundary,¶ms4d) + IceRayTracing::fpathD(-(IceRayTracing::TransitionBoundary+1e-7),¶ms4f) - IceRayTracing::fpathD(z0,¶ms4a); + pathRa2[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(z1,¶ms4b); + } + }else{ + timeRa1[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(z0,¶ms4a); + timeRa2[i]=IceRayTracing::ftimeD(-zmax[i],¶ms4c) - IceRayTracing::ftimeD(z1,¶ms4b); + + pathRa1[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(z0,¶ms4a); + pathRa2[i]=IceRayTracing::fpathD(-zmax[i],¶ms4c) - IceRayTracing::fpathD(z1,¶ms4b); + } + + raytime[i]= timeRa1[i] + timeRa2[i]; + pathRa[i]= pathRa1[i] + pathRa2[i]; + + if(Flip==true){ + double dumRa=timeRa2[i]; + timeRa2[i]=timeRa1[i]; + timeRa1[i]=dumRa; + + dumRa=pathRa2[i]; + pathRa2[i]=pathRa1[i]; + pathRa1[i]=dumRa; + } + } + timeRa[i]=raytime[i]; + + /* Setup the function that will be used to calculate the angle of reception for all the rays */ + gsl_function F5; + struct IceRayTracing::fDnfR_params params5c = {IceRayTracing::A_ice, IceRayTracing::GetB(z1), IceRayTracing::GetC(z1), lvalueRa[i]}; + double result, abserr; + F5.function = &IceRayTracing::fDnfR; + + /* Calculate the recieve angle for refacted ray by calculating the derivative of the function at the Rx position */ + F5.params = ¶ms5c; + gsl_deriv_central (&F5, z1, 1e-8, &result, &abserr); + RangRa[i]=180-atan(result)*(180.0/IceRayTracing::pi); + + /* When the Tx and Rx are at the same depth my function struggles to find a ray between them when they are very close to each other. In that case the ray is pretty much like a straight line. */ + if(z1==z0 && std::isnan(RangRa[i])==true){ + RangRa[i]=180-LangRa[i]; + } + + /* This sometimes happens that when the Rx is very close to the peak point (or the turning point) of the ray then its hard to calculate the derivative around that area since the solution blows up around that area. therefore this is a good approximation. */ + if(z1!=z0 && std::isnan(RangRa[i])==true){ + RangRa[i]=90; + } + + }//// i loop + + dsw=0; + /* If the Tx and Rx depth were switched then put them back to their original position */ + if(Flip==true){ + dsw=z0; + z0=z1; + z1=dsw; + } + + output[0]=RangRa[0]; + output[1]=LangRa[0]; + output[2]=timeRa[0]; + output[3]=lvalueRa[0]; + output[4]=checkzeroRa[0]; + output[5]=timeRa1[0]; + output[6]=timeRa2[0]; + output[7]=zmax[0]; + + /* If the flip case is true where we flipped Rx and Tx depths to trace rays then make sure everything is switched back before we give the output to the user. */ + if(Flip==true){ + output[0]=180-LangRa[0]; + output[1]=180-RangRa[0]; + } + + output[8]=RangRa[1]; + output[9]=LangRa[1]; + output[10]=timeRa[1]; + output[11]=lvalueRa[1]; + output[12]=checkzeroRa[1]; + output[13]=timeRa1[1]; + output[14]=timeRa2[1]; + output[15]=zmax[1]; + + output[16]=pathRa[0]; + output[17]=pathRa1[0]; + output[18]=pathRa2[0]; + + output[19]=pathRa[1]; + output[20]=pathRa1[1]; + output[21]=pathRa2[1]; + + /* If the flip case is true where we flipped Rx and Tx depths to trace rays then make sure everything is switched back before we give the output to the user. */ + if(Flip==true){ + output[8]=180-LangRa[1]; + output[9]=180-RangRa[1]; + } + + return output; +} + + +/* This function returns the x and z values for the full Direct ray and prints out the ray path in a text file */ +void IceRayTracing::GetFullDirectRayPath(double z0, double x1, double z1,double lvalueD, vector &x, vector &z){ + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* Set the name of the text files */ + //ofstream aoutD("DirectRay.txt"); + + /* Set the step size for plotting */ + double h=0.5; + /* Set the total steps required for looping over the whole ray path */ + int dmax=100000; + /* Set the values to start the rays from */ + double zn=z1; + double xn=0; + + /* Map out the direct ray path */ + int npnt=0; + double checknan=0; + struct IceRayTracing::fDnfR_params params6a; + struct IceRayTracing::fDnfR_params params6b; + struct IceRayTracing::fDnfR_params params6c; + struct IceRayTracing::fDnfR_params params6d; + + + for(int i=0;iIceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)>IceRayTracing::TransitionBoundary){ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(z0,¶ms6b); + } + if (fabs(z0)==IceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)==IceRayTracing::TransitionBoundary){ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(-IceRayTracing::TransitionBoundary,¶ms6c) + IceRayTracing::fDnfR(-(IceRayTracing::TransitionBoundary+1e-7),¶ms6d) - IceRayTracing::fDnfR(z0,¶ms6b); + } + }else{ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(z0,¶ms6b); + } + + checknan=IceRayTracing::fDnfR(zn,¶ms6a); + if(std::isnan(checknan)==false && Flip==false){ + x.push_back(xn); + z.push_back(zn); + //aoutD< &x, vector &z){ + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* Set the name of the text files */ + //ofstream aoutR("ReflectedRay.txt"); + /* Set the step size for plotting. */ + double h=0.5; + /* Set the total steps required for looping over the whole ray path */ + int dmax=100000; + /* Set the values to start the rays from */ + double zn=z1; + double xn=0; + + /* Map out the direct ray path */ + int npnt=0; + double checknan=0; + struct IceRayTracing::fDnfR_params params6a; + struct IceRayTracing::fDnfR_params params6b; + struct IceRayTracing::fDnfR_params params6c; + struct IceRayTracing::fDnfR_params params6f; + struct IceRayTracing::fDnfR_params params6d; + + /* Map out the 1st part of the reflected ray */ + for(int i=0;iIceRayTracing::TransitionBoundary && fabs(zn)>IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(1e-7,¶ms6c) - IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary,¶ms6d) + IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary+1e-7,¶ms6f) - IceRayTracing::fDnfR(-z0,¶ms6b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)==IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary,¶ms6d) + IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary+1e-7,¶ms6f) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(1e-7,¶ms6c) - IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary,¶ms6d) + IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary+1e-7,¶ms6f) - IceRayTracing::fDnfR(-z0,¶ms6b); + } + }else{ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(1e-7,¶ms6c) - IceRayTracing::fDnfR(-z0,¶ms6b); + } + + xn= distancez0z1 - 2*(distancez0surface); + + checknan=IceRayTracing::fDnfR(-zn,¶ms6a); + if(std::isnan(checknan)==false && zn<=0 && Flip==false){ + x.push_back(xn); + z.push_back(zn); + //aoutR<0){ + i=dmax+2; + } + } + + /* Map out the 2nd part of the reflected ray */ + zn=-1e-7; + for(int i=0;iIceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)>IceRayTracing::TransitionBoundary){ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(z0,¶ms6b); + } + if (fabs(z0)==IceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)==IceRayTracing::TransitionBoundary){ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(-IceRayTracing::TransitionBoundary,¶ms6c) + IceRayTracing::fDnfR(-(IceRayTracing::TransitionBoundary+1e-6),¶ms6d) - IceRayTracing::fDnfR(z0,¶ms6b); + } + }else{ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(z0,¶ms6b); + } + + checknan=IceRayTracing::fDnfR(zn,¶ms6a); + if(std::isnan(checknan)==false && Flip==false){ + x.push_back(xn); + z.push_back(zn); + //aoutR< &x, vector &z,int raynumber){ + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* Set the name of the text files */ + //ofstream aoutRa; + /* Set the name of the text files */ + // if(raynumber==1){ + // aoutRa.open("RefractedRay1.txt"); + // } + // if(raynumber==2){ + // aoutRa.open("RefractedRay2.txt"); + // } + /* Set the step size for plotting. */ + double h=0.5; + /* Set the total steps required for looping over the whole ray path */ + int dmax=100000; + /* Set the values to start the rays from */ + double zn=z1; + double xn=0; + + /* Map out the direct ray path */ + int npnt=0; + double checknan=0; + struct IceRayTracing::fDnfR_params params6a; + struct IceRayTracing::fDnfR_params params6b; + struct IceRayTracing::fDnfR_params params6c; + struct IceRayTracing::fDnfR_params params6f; + struct IceRayTracing::fDnfR_params params6d; + + /* Map out the 1st part of the refracted ray */ + for(int i=0;iIceRayTracing::TransitionBoundary && fabs(zn)>IceRayTracing::TransitionBoundary){ + if(zmax<=IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(zmax,¶ms6c) - IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary,¶ms6d) + IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary+1e-6,¶ms6f) - IceRayTracing::fDnfR(-z0,¶ms6b); + }else{ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(zmax,¶ms6c) - IceRayTracing::fDnfR(-z0,¶ms6b); + } + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)==IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary,¶ms6d) + IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary+1e-6,¶ms6f) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(zmax,¶ms6c) - IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary,¶ms6d) + IceRayTracing::fDnfR(IceRayTracing::TransitionBoundary+1e-6,¶ms6f) - IceRayTracing::fDnfR(-z0,¶ms6b); + } + }else{ + distancez0z1=IceRayTracing::fDnfR(-zn,¶ms6a) - IceRayTracing::fDnfR(-z0,¶ms6b); + distancez0surface=IceRayTracing::fDnfR(zmax,¶ms6c) - IceRayTracing::fDnfR(-z0,¶ms6b); + } + + xn= distancez0z1 - 2*(distancez0surface); + + checknan=IceRayTracing::fDnfR(-zn,¶ms6a); + if(std::isnan(checknan)==false && zn<=0 && Flip==false){ + x.push_back(xn); + z.push_back(zn); + //aoutRa<-zmax){ + i=dmax+2; + } + } + + /* Map out the 2nd part of the refracted ray */ + zn=-zmax; + for(int i=0;iIceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)>IceRayTracing::TransitionBoundary){ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(z0,¶ms6b); + } + if (fabs(z0)==IceRayTracing::TransitionBoundary && fabs(zn)IceRayTracing::TransitionBoundary && fabs(zn)==IceRayTracing::TransitionBoundary){ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(-IceRayTracing::TransitionBoundary,¶ms6c) + IceRayTracing::fDnfR(-(IceRayTracing::TransitionBoundary+1e-6),¶ms6d) - IceRayTracing::fDnfR(z0,¶ms6b); + } + }else{ + xn=IceRayTracing::fDnfR(zn,¶ms6a) - IceRayTracing::fDnfR(z0,¶ms6b); + } + + checknan=IceRayTracing::fDnfR(zn,¶ms6a); + if(std::isnan(checknan)==false && Flip==false){ + x.push_back(xn); + z.push_back(zn); + //aoutRa< xD,zD; + vector xR,zR; + vector xRa[2],zRa[2]; + GetFullDirectRayPath(z0,x1,z1,lvalueD,xD,zD); + GetFullReflectedRayPath(z0,x1,z1,lvalueR,xR,zR); + if((fabs(checkzeroR)>0.5 || fabs(checkzeroD)>0.5) && fabs(checkzeroRa[0])<0.5){ + GetFullRefractedRayPath(z0,x1,z1,zmax[0],lvalueRa[0],xRa[0],zRa[0],1); + if(fabs(checkzeroRa[1])<0.5){ + GetFullRefractedRayPath(z0,x1,z1,zmax[1],lvalueRa[1],xRa[1],zRa[1],2); + } + } + + /* An example of how the function returns the filled in vector and you can see the printed out values as a sanity check */ + // for (int i=0;i0.5 || fabs(checkzeroD)>0.5){ + double* GetRefractedRay=IceRayTracing::GetRefractedRayPar(z0,x1,z1,LangR,RangR,checkzeroD,checkzeroR); + RangRa[0]=GetRefractedRay[0]; + LangRa[0]=GetRefractedRay[1]; + timeRa[0]=GetRefractedRay[2]; + lvalueRa[0]=GetRefractedRay[3]; + checkzeroRa[0]=GetRefractedRay[4]; + timeRa1[0]=GetRefractedRay[5]; + timeRa2[0]=GetRefractedRay[6]; + zmax[0]=GetRefractedRay[7]; + + if(fabs(checkzeroR)>0.5 && fabs(checkzeroD)>0.5){ + RangRa[1]=GetRefractedRay[8]; + LangRa[1]=GetRefractedRay[9]; + timeRa[1]=GetRefractedRay[10]; + lvalueRa[1]=GetRefractedRay[11]; + checkzeroRa[1]=GetRefractedRay[12]; + timeRa1[1]=GetRefractedRay[13]; + timeRa2[1]=GetRefractedRay[14]; + zmax[1]=GetRefractedRay[15]; + } + + pathRa[0]=GetRefractedRay[16]; + pathRa1[0]=GetRefractedRay[17]; + pathRa2[0]=GetRefractedRay[18]; + + pathRa[1]=GetRefractedRay[19]; + pathRa1[1]=GetRefractedRay[20]; + pathRa2[1]=GetRefractedRay[21]; + + delete []GetRefractedRay; + } + + /* This part of the code can be used if the user wants to plot the individual ray paths. This part of the code prints out the individual ray paths in text files and also plots them on a canvas */ + if(PlotRayPaths==true){ + double lvalues[4]; + lvalues[0]=lvalueD; + lvalues[1]=lvalueR; + lvalues[2]=lvalueRa[0]; + lvalues[3]=lvalueRa[1]; + + double checkzeroes[4]; + checkzeroes[0]=checkzeroD; + checkzeroes[1]=checkzeroR; + checkzeroes[2]=checkzeroRa[0]; + checkzeroes[3]=checkzeroRa[1]; + + IceRayTracing::PlotAndStoreRays(x0,z0,z1,x1,zmax,lvalues,checkzeroes); + } + + /* print out all the output from the code */ + //cout<<0<<" ,x0= "<0.5){ + output[8]=-1000; + } + if(fabs(checkzeroR)>0.5){ + output[9]=-1000; + } + if(fabs(checkzeroRa[0])>0.5){ + output[10]=-1000; + } + if(fabs(checkzeroRa[1])>0.5){ + output[11]=-1000; + } + + return output; +} + +/* Analytical solution describing ray paths in ice as function of depth for constant refractive index*/ +double IceRayTracing::fDnfR_Cnz(double x,void *params){ + + struct IceRayTracing::fDnfR_params *p= (struct IceRayTracing::fDnfR_params *) params; + double A = p->a; + double L = p->l; + + return (L/sqrt(A*A-L*L))*x; +} + +/* Analytical solution describing the ray path in ice as a function of the L parameter for constant refractive index*/ +double IceRayTracing::fDnfR_L_Cnz(double x,void *params){ + + struct IceRayTracing::fDnfR_L_params *p= (struct IceRayTracing::fDnfR_L_params *) params; + double A = p->a; + double Z = p->z; + + double out=0; + if(A>x){ + out=(x/sqrt(A*A-x*x))*Z; + }else{ + out=tan(asin(x/A))*Z; + } + return out; +} + +/* This function is minimised to find the launch angle (or the L parameter) for the reflected ray for constant refractive index*/ +double IceRayTracing::fRa_Cnz(double x,void *params){ + struct IceRayTracing::fDanfRa_params *p= (struct IceRayTracing::fDanfRa_params *) params; + double A = p->a; + double z0 = p->z0; + double x1 = p->x1; + double z1 = p->z1; + + struct IceRayTracing::fDnfR_L_params params1a = {A, 0, 0, -z1}; + struct IceRayTracing::fDnfR_L_params params1b = {A, 0, 0, -z0}; + struct IceRayTracing::fDnfR_L_params params1c = {A, 0, 0, 0.0}; + + return IceRayTracing::fDnfR_L_Cnz(x,¶ms1a) - IceRayTracing::fDnfR_L_Cnz(x,¶ms1b) - 2*( IceRayTracing::fDnfR_L_Cnz(x,¶ms1c) - IceRayTracing::fDnfR_L_Cnz(x,¶ms1b) ) - x1; +} + +double IceRayTracing::fRa_Cnz_df(double x,void *params){ + gsl_function F; + F.function = &IceRayTracing::fRa_Cnz; + F.params = params; + + double result,abserr; + gsl_deriv_central (&F, x, 1e-8, &result, &abserr); + + return result; +} + +void IceRayTracing::fRa_Cnz_fdf (double x, void *params,double *y, double *dy){ + *y = IceRayTracing::fRa_Cnz(x,params); + *dy = IceRayTracing::fRa_Cnz_df(x,params); +} + +/* This functions works for the Direct ray and gives you back the launch angle, receive angle and propagation time of the ray together with values of the L parameter. This for constant refractive index*/ +double* IceRayTracing::GetDirectRayPar_Cnz(double z0, double x1, double z1, double A_ice_Cnz){ + + double *output=new double[4]; + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* Calculate the launch angle and the value of the L parameter */ + double LangD=(IceRayTracing::pi*0.5-atan(fabs(z1-z0)/x1))*(180.0/IceRayTracing::pi); + double lvalueD=A_ice_Cnz*sin(LangD*(IceRayTracing::pi/180.0)); + double timeD=(sqrt( pow(x1,2) + pow(z1-z0,2) )/IceRayTracing::c_light_ms)*A_ice_Cnz; + + /* Calculate the recieve angle for direct rays by which is the same as the launch angle */ + double RangD=LangD; + + dsw=0; + /* If the Tx and Rx depth were switched then put them back to their original position */ + if(Flip==true){ + dsw=z0; + z0=z1; + z1=dsw; + } + + output[0]=RangD; + output[1]=LangD; + output[2]=timeD; + output[3]=lvalueD; + + /* If the flip case is true where we flipped Rx and Tx depths to trace rays then make sure everything is switched back before we give the output to the user. */ + if(Flip==true){ + output[0]=180-LangD; + output[1]=180-RangD; + } + + return output; +} + +/* This functions works for the Reflected ray and gives you back the launch angle, receive angle and propagation times (of the whole ray and the two direct rays that make it up) together with values of the L parameter. This is for constant refractive index*/ +double *IceRayTracing::GetReflectedRayPar_Cnz(double z0, double x1 , double z1, double A_ice_Cnz){ + + double *output=new double[8]; + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* First we setup the fRa function that will be minimised to get the launch angle (or the L parameter) for the reflected ray. */ + gsl_function F3; + struct IceRayTracing::fDanfRa_params params3= {A_ice_Cnz, z0, x1, z1}; + F3.function = &IceRayTracing::fRa_Cnz; + F3.params = ¶ms3; + + // gsl_function_fdf F3; + // struct IceRayTracing::fDanfRa_params params3= {A_ice_Cnz, z0, x1, z1}; + // F3.f = &fRa_Cnz; + // F3.df = &fRa_Cnz_df; + // F3.fdf = &fRa_Cnz_fdf; + // F3.params = ¶ms3; + + /* In my raytracing solution given in the function IceRayTracing::fDnfR_Cnz the launch angle (or the L parameter) has limit placed on it by this part in the solution that L &x, vector &z){ + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + /* Set the name of the text files */ + //ofstream aoutD("DirectRay_Cnz.txt"); + /* Set the step size for plotting */ + double h=0.5; + /* Set the total steps required for looping over the whole ray path */ + int dmax=100000; + /* Set the values to start the rays from */ + double zn=z1; + double xn=0; + + /* Map out the direct ray path */ + int npnt=0; + double checknan=0; + struct IceRayTracing::fDnfR_params params6a; + struct IceRayTracing::fDnfR_params params6b; + + for(int i=0;i &x, vector &z){ + + /* My raytracer can only work the Tx is below the Rx. If the Tx is higher than the Rx than we need to flip the depths to allow for raytracing and then we will flip them back later at the end */ + bool Flip=false; + double dsw=z0; + if(z0>z1){ + z0=z1; + z1=dsw; + Flip=true; + } + + // /* Set the name of the text files */ + //ofstream aoutR("ReflectedRay_Cnz.txt"); + /* Set the step size for plotting. */ + double h=0.5; + /* Set the total steps required for looping over the whole ray path */ + int dmax=100000; + /* Set the values to start the rays from */ + double zn=z1; + double xn=0; + + /* Map out the direct ray path */ + int npnt=0; + double checknan=0; + struct IceRayTracing::fDnfR_params params6a; + struct IceRayTracing::fDnfR_params params6b; + struct IceRayTracing::fDnfR_params params6c; + + /* Map out the 1st part of the reflected ray */ + for(int i=0;i0){ + i=dmax+2; + } + } + + /* Map out the 2nd part of the reflected ray */ + zn=0.0; + for(int i=0;i xD,zD; + vector xR,zR; + GetFullDirectRayPath_Cnz(z0,x1,z1,lvalueD,A_ice_Cnz,xD,zD); + GetFullReflectedRayPath_Cnz(z0,x1,z1,lvalueR,A_ice_Cnz,xR,zR); + +} + +/* This is the main raytracing function. x0 always has to be zero. z0 is the Tx depth in m and z1 is the depth of the Rx in m. Both depths are negative. x1 is the distance between them. This functions works for a constant refractive index */ +double *IceRayTracing::IceRayTracing_Cnz(double x0, double z0, double x1, double z1, double A_ice_Cnz){ + + /* define a pointer to give back the output of raytracing */ + double *output=new double[9]; + + /* Store the ray paths in text files */ + bool PlotRayPaths=false; + + /* ********This part of the code will try to get the Direct ray between Rx and Tx.********** */ + double* GetDirectRay=GetDirectRayPar_Cnz(z0,x1,z1,A_ice_Cnz); + double RangD=GetDirectRay[0]; + double LangD=GetDirectRay[1]; + double timeD=GetDirectRay[2]; + double lvalueD=GetDirectRay[3]; + double checkzeroD=GetDirectRay[4]; + delete []GetDirectRay; + + /* ********This part of the code will try to get the Reflected ray between Rx and Tx.********** */ + double* GetReflectedRay=GetReflectedRayPar_Cnz(z0,x1,z1,A_ice_Cnz); + double RangR=GetReflectedRay[0]; + double LangR=GetReflectedRay[1]; + double timeR=GetReflectedRay[2]; + double lvalueR=GetReflectedRay[3]; + double checkzeroR=GetReflectedRay[4]; + double timeR1=GetReflectedRay[5]; + double timeR2=GetReflectedRay[6]; + double AngleOfIncidenceInIce=GetReflectedRay[7]; + delete []GetReflectedRay; + + /* This part of the code can be used if the user wants to plot the individual ray paths. This part of the code prints out the individual ray paths in text files and also plots them on a canvas */ + if(PlotRayPaths==true){ + double lvalues[2]; + lvalues[0]=lvalueD; + lvalues[1]=lvalueR; + + PlotAndStoreRays_Cnz(x0,z0,z1,x1,lvalues,A_ice_Cnz); + } + + /* Fill in the output pointer after calculating all the results */ + output[0]=LangD; + output[1]=LangR; + output[2]=timeD; + output[3]=timeR; + output[4]=RangD; + output[5]=RangR; + + /* fill in the output array part where you fill in the times for the two parts of the reflected or refracted rays */ + output[6]=timeR1; + output[7]=timeR2; + output[8]=AngleOfIncidenceInIce; + + + return output; +} + +/* The set of functions starting with the name "fDa" are used in the minimisation procedure to find the launch angle (or the L parameter) for the direct ray */ +double IceRayTracing::fDa_Air(double x,void *params){ + struct IceRayTracing::fDanfRa_params *p= (struct IceRayTracing::fDanfRa_params *) params; + double A = p->a; + double z0 = p->z0; + double x1 = p->x1; + double z1 = p->z1; + + double nz_Air=1; + double AngleInAir=asin(x/nz_Air); + double x1_Air=z1*tan(AngleInAir); + + struct IceRayTracing::fDnfR_L_params params1a = {A, IceRayTracing::GetB(1e-7), IceRayTracing::GetC(1e-7), -1e-7}; + struct IceRayTracing::fDnfR_L_params params1b = {A, IceRayTracing::GetB(z0), IceRayTracing::GetC(z0), z0}; + struct IceRayTracing::fDnfR_L_params params1c = {A, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), IceRayTracing::GetC(IceRayTracing::TransitionBoundary), -IceRayTracing::TransitionBoundary}; + struct IceRayTracing::fDnfR_L_params params1d = {A, IceRayTracing::GetB(-(IceRayTracing::TransitionBoundary+0.000001)), IceRayTracing::GetC(-(IceRayTracing::TransitionBoundary+0.000001)), -(IceRayTracing::TransitionBoundary+0.000001)}; + + double distancez0z1=0; + + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(1e-7)>IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(1e-7)IceRayTracing::TransitionBoundary && fabs(1e-7)==IceRayTracing::TransitionBoundary){ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1c) + IceRayTracing::fDnfR_L(x,¶ms1d) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + }else{ + distancez0z1=IceRayTracing::fDnfR_L(x,¶ms1a) - IceRayTracing::fDnfR_L(x,¶ms1b); + } + // if(std::isnan(distancez0z1)){ + // distancez0z1=1e9; + // } + if(std::isnan(x1_Air)){ + x1_Air=1e9; + } + + double output=distancez0z1+x1_Air-x1; + + return output; +} + +/* This functions works for the Direct ray and gives you back the launch angle, receive angle and propagation time of the ray together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. */ +double* IceRayTracing::GetDirectRayPar_Air(double z0, double x1, double z1){ + + double *output=new double[5]; + + /* First we setup the fDa function that will be minimised to get the launch angle (or the L parameter) for the direct ray. */ + gsl_function F1; + struct IceRayTracing::fDanfRa_params params1= {IceRayTracing::A_ice, z0, x1, z1}; + F1.function = &IceRayTracing::fDa_Air; + F1.params = ¶ms1; + + /* In my raytracing solution given in the function IceRayTracing::fDnfR the launch angle (or the L parameter) has limit placed on it by this part in the solution sqrt( n(z)^2 - L^2) . This sqrt cannot be negative for both z0 and 1e-7 and this sets the upper limit in our minimisation to get the launch angle (or the L parameter). Here I am basically setting the upper limit as GSL requires that my function is well behaved on the upper and lower bounds I give it for minimisation. */ + double UpLimnz[]={IceRayTracing::Getnz(1e-7),IceRayTracing::Getnz(z0)}; + double* UpperLimitL=min_element(UpLimnz,UpLimnz+2); + + /* Do the minimisation and get the value of the L parameter and the launch angle and then verify to see that the value of L that we got was actually a root of fDa function. */ + double lvalueD=IceRayTracing::FindFunctionRoot(F1,1e-7,UpperLimitL[0]); + double LangD=asin(lvalueD/IceRayTracing::Getnz(z0))*(180.0/IceRayTracing::pi); + double checkzeroD=IceRayTracing::fDa_Air(lvalueD,¶ms1); + + /* Get the propagation time for the direct ray using the ftimeD function after we have gotten the value of the L parameter. */ + struct IceRayTracing::ftimeD_params params2a = {IceRayTracing::A_ice, IceRayTracing::GetB(z0), -IceRayTracing::GetC(z0), IceRayTracing::c_light_ms,lvalueD}; + struct IceRayTracing::ftimeD_params params2b = {IceRayTracing::A_ice, IceRayTracing::GetB(1e-7), -IceRayTracing::GetC(1e-7), IceRayTracing::c_light_ms,lvalueD}; + struct IceRayTracing::ftimeD_params params2c = {IceRayTracing::A_ice, IceRayTracing::GetB(IceRayTracing::TransitionBoundary), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary), IceRayTracing::c_light_ms, lvalueD}; + struct IceRayTracing::ftimeD_params params2d = {IceRayTracing::A_ice, IceRayTracing::GetB(IceRayTracing::TransitionBoundary+1e-7), -IceRayTracing::GetC(IceRayTracing::TransitionBoundary+1e-7), IceRayTracing::c_light_ms, lvalueD}; + + /* we do the subtraction because we are measuring the time taken between the Tx and Rx positions */ + double timeD=0; + if(IceRayTracing::TransitionBoundary!=0){ + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(1e-7)>IceRayTracing::TransitionBoundary){ + timeD=IceRayTracing::ftimeD(-z0,¶ms2a) - IceRayTracing::ftimeD(1e-7,¶ms2b); + } + if (fabs(z0)>IceRayTracing::TransitionBoundary && fabs(1e-7)IceRayTracing::TransitionBoundary && fabs(1e-7)==IceRayTracing::TransitionBoundary){ + timeD=IceRayTracing::ftimeD(-z0,¶ms2a) - IceRayTracing::ftimeD(IceRayTracing::TransitionBoundary+1e-7,¶ms2d) + IceRayTracing::ftimeD(IceRayTracing::TransitionBoundary,¶ms2c) - IceRayTracing::ftimeD(1e-7,¶ms2b); + } + }else{ + timeD=IceRayTracing::ftimeD(-z0,¶ms2a) - IceRayTracing::ftimeD(1e-7,¶ms2b); + } + + /* Setup the function that will be used to calculate the angle of reception for all the rays */ + gsl_function F5; + struct IceRayTracing::fDnfR_params params5a = {IceRayTracing::A_ice, IceRayTracing::GetB(1e-7), -IceRayTracing::GetC(1e-7), lvalueD}; + double result, abserr; + F5.function = &IceRayTracing::fDnfR; + + /* Calculate the recieve angle for direc rays by calculating the derivative of the function at the Rx position */ + F5.params = ¶ms5a; + gsl_deriv_central (&F5, 1e-7, 1e-8, &result, &abserr); + double RangD=atan(result); + + /* When the Tx and Rx are at the same depth my function struggles to find a ray between them when they are very close to each other. In that case the ray is pretty much like a straight line. */ + if(z1==z0 && std::isnan(RangD)==true){ + RangD=180-LangD; + } + + /* This sometimes happens that when the Rx is very close to the peak point (or the turning point) of the ray then its hard to calculate the derivative around that area since the solution blows up around that area. therefore this is a good approximation. */ + if(z1!=z0 && std::isnan(RangD)==true){ + RangD=90; + } + + double AirAngle=asin(IceRayTracing::Getnz(1e-7)*sin(RangD)); + double AirHorizontalDistance=tan(AirAngle)*z1; + double AirTime=AirHorizontalDistance/IceRayTracing::c_light_ms; + + timeD=timeD+AirTime; + RangD=AirAngle*(180.0/IceRayTracing::pi); + + output[0]=RangD; + output[1]=LangD; + output[2]=timeD; + output[3]=lvalueD; + output[4]=checkzeroD; + + if(fabs(checkzeroD)>0.5){ + output[0]=-1000; + } + + return output; +} + +double *IceRayTracing::DirectRayTracer(double xT, double yT, double zT, double xR, double yR, double zR){ + double TxCor[3]={xT,yT,zT}; + double RxCor[3]={xR,yR,zR}; + + ////For recording how much time the process took + //auto t1b = std::chrono::high_resolution_clock::now(); + + double x0=0;/////always has to be zero + double z0=TxCor[2]; + double x1=sqrt(pow(TxCor[0]-RxCor[0],2)+pow(TxCor[1]-RxCor[1],2)); + double z1=RxCor[2]; + + double * getresults=IceRayTracing::IceRayTracing(x0,z0,x1,z1); + + // cout<<"*******For the Direct Ray********"< OutputValues[5]; + + if(getresults[8]!=-1000){ + // cout<<"*******For the Direct Ray********"<( t2b - t1b ).count(); + + // double Duration=durationb; + //cout<<"total time taken by the script: "<0){ + IceRayTracing::GridStartZ=-20; + IceRayTracing::GridStopZ=0; + //} + + //cout<<"Grid Variables are "< xRay[2]; + vector zRay[2]; + IceRayTracing::GetRayTracingSolutions(zR, xR, zT, TimeRay_Tx, PathRay_Tx, LaunchAngle_Tx, RecieveAngle_Tx, IgnoreCh_Tx, IncidenceAngleInIce_Tx, A0, frequency, AttRay_Tx,xRay,zRay); + + GridPositionXb[AntNum][ix]=xR; + GridPositionZb[AntNum][iz]=zR; + + if(IgnoreCh_Tx[0]!=0){ + GridZValueb[AntNum][0].push_back(TimeRay_Tx[0]); + GridZValueb[AntNum][1].push_back(PathRay_Tx[0]); + GridZValueb[AntNum][2].push_back(LaunchAngle_Tx[0]); + GridZValueb[AntNum][3].push_back(RecieveAngle_Tx[0]); + GridZValueb[AntNum][4].push_back(AttRay_Tx[0]); + }else{ + GridZValueb[AntNum][0].push_back(-1000); + GridZValueb[AntNum][1].push_back(-1000); + GridZValueb[AntNum][2].push_back(-1000); + GridZValueb[AntNum][3].push_back(-1000); + GridZValueb[AntNum][4].push_back(-1000); + } + + if(IgnoreCh_Tx[1]!=0){ + GridZValueb[AntNum][5].push_back(TimeRay_Tx[1]); + GridZValueb[AntNum][6].push_back(PathRay_Tx[1]); + GridZValueb[AntNum][7].push_back(LaunchAngle_Tx[1]); + GridZValueb[AntNum][8].push_back(RecieveAngle_Tx[1]); + GridZValueb[AntNum][9].push_back(AttRay_Tx[1]); + }else{ + GridZValueb[AntNum][5].push_back(-1000); + GridZValueb[AntNum][6].push_back(-1000); + GridZValueb[AntNum][7].push_back(-1000); + GridZValueb[AntNum][8].push_back(-1000); + GridZValueb[AntNum][9].push_back(-1000); + } + + } + } + + auto t2b = std::chrono::high_resolution_clock::now(); + double duration = std::chrono::duration_cast( t2b - t1b ).count(); + cout<<"The table took "<IceRayTracing::TotalStepsX_O){ + minXbin=IceRayTracing::TotalStepsX_O-2; + } + if(minZbin+1>IceRayTracing::TotalStepsZ_O){ + minZbin=IceRayTracing::TotalStepsZ_O-2; + } + + int startbinX=minXbin-1; + int endbinX=minXbin+1; + int startbinZ=minZbin-1; + int endbinZ=minZbin+1; + + newXbin=(minXbin-1); + newZbin=(minZbin-1); + //int newich=newXbin*IceRayTracing::TotalStepsZ_O+newZbin; + double minDist1=fabs(((xR-GridPositionXb[AntNum][newXbin])*(xR-GridPositionXb[AntNum][newXbin])+(zR-GridPositionZb[AntNum][newZbin])*(zR-GridPositionZb[AntNum][newZbin]))); + + newXbin=minXbin+1; + newZbin=minZbin+1; + //newich=newXbin*IceRayTracing::TotalStepsZ_O+newZbin; + double minDist2=fabs(((xR-GridPositionXb[AntNum][newXbin])*(xR-GridPositionXb[AntNum][newXbin])+(zR-GridPositionZb[AntNum][newZbin])*(zR-GridPositionZb[AntNum][newZbin]))); + + startbinX=minXbin-1; + endbinX=minXbin+1; + startbinZ=minZbin-1; + endbinZ=minZbin+1; + + sum1=0; + sum2=0; + NewZValue=0; + + for(int ixn=startbinX;ixn=0 && newich=0 && izn>=0){ + MinDist[count]=fabs(((xR-GridPositionXb[AntNum][ixn])*(xR-GridPositionXb[AntNum][ixn])+(zR-GridPositionZb[AntNum][izn])*(zR-GridPositionZb[AntNum][izn]))); + MinDistBin[count]=newich; + + if(GridZValueb[AntNum][rtParameter][MinDistBin[count]]!=-1000){ + //cout<<"in here "< xRay[2], vector zRay[2]){ + + int DHits=0,RHits=0; + double timeD, timeR, timeRa[2]; + double pathD, pathR, pathRa[2]; + double RangD, RangR, RangRa[2]; + double LangD, LangR, LangRa[2]; + double LvalueD, LvalueR, LvalueRa[2],zmax[2]; + double AttD, AttR, AttRa[2]; + int RayType[2]; + + // double Distance=sqrt(pow(RxCor[0]-TxCor[0],2) + pow(RxCor[1]-TxCor[1],2)); + // double RxDepth=AntennaCoordRx[2]+AvgAntennaCoordRx[2]; + // double TxDepth=AntennaCoordTx[2]; + //cout<<"rt parameters are "<<0<<" "<TimeRay[1] && RecieveAngle[0]!=-1000 && RecieveAngle[1]!=-1000){ + swap(LaunchAngle[0],LaunchAngle[1]); + swap(RecieveAngle[0],RecieveAngle[1]); + swap(TimeRay[0],TimeRay[1]); + swap(RayType[0],RayType[1]); + swap(PathRay[0],PathRay[1]); + swap(AttRay[0],AttRay[1]); + } + + //vector xRay[2]; + //vector zRay[2]; + + for(int iray=0;iray<2;iray++){ + if(IgnoreCh[iray]!=0){ + if(RayType[iray]==1 ){ + IceRayTracing::GetFullDirectRayPath(TxDepth,Distance,RxDepth,LvalueD,xRay[iray],zRay[iray]); + } + if(RayType[iray]==2){ + IceRayTracing::GetFullReflectedRayPath(TxDepth,Distance,RxDepth,LvalueR,xRay[iray],zRay[iray]); + } + if(RayType[iray]==3){ + IceRayTracing::GetFullRefractedRayPath(TxDepth,Distance,RxDepth,zmax[0],LvalueRa[0],xRay[iray],zRay[iray],1); + } + if(RayType[iray]==4){ + IceRayTracing::GetFullRefractedRayPath(TxDepth,Distance,RxDepth,zmax[1],LvalueRa[1],xRay[iray],zRay[iray],2); + } + } + } + + // cout<<"inside Ray 1 "< +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace std; + +namespace IceRayTracing{ + + /********Stuff for Interpolation**********/ + static vector> GridPositionXb; + static vector> GridPositionZb; + static vector>> GridZValueb; + + static double GridStepSizeX_O=0.1; + static double GridStepSizeZ_O=0.1; + static double GridWidthX=20; + static double GridWidthZ=20; + + static int GridPoints=100;////just set a non-zeronumber for now + static int TotalStepsX_O=100;////just set a non-zeronumber for now + static int TotalStepsZ_O=100;////just set a non-zeronumber for now + static double GridStartX=40;////just set a non-zeronumber for now + static double GridStopX=60;////just set a non-zeronumber for now + static double GridStartZ=-20;////just set a non-zeronumber for now + static double GridStopZ=0; + + /* Set the value of pi */ + static constexpr double pi=3.14159265359; + /* Set the value of the speed of light in m/s */ + static constexpr double c_light_ms=299792458; + /* Set the value of the asymptotic parameter of the refractive index model */ + + static const double A_ice_def=1.78; + static const double B_ice_def=-0.43; + static const double C_ice_def=0.0132; + + static double A_ice=A_ice_def; + static double B_ice=B_ice_def; + static double C_ice=C_ice_def; + static constexpr double TransitionBoundary=0; + // const double A_ice=1.775; + // const double TransitionBoundary=14.9; + + + /* Get the value of the B parameter for the refractive index model */ + void SetA(double &A); + + /* Get the value of the B parameter for the refractive index model */ + void SetB(double &B); + + /* Get the value of the C parameter for the refractive index model */ + void SetC(double &C); + + /* Get the value of the B parameter for the refractive index model */ + double GetB(double z); + + /* Get the value of the C parameter for the refractive index model */ + double GetC(double z); + + /* Get the value of refractive index model for a given depth */ + double Getnz(double z); + + /* E-feild Power Fresnel coefficient for S-polarised wave which is perpendicular to the plane of propogation/incidence. This function gives you back the reflectance. The transmittance is T=1-R */ + double Refl_S(double thetai); + +/* E-feild Power Fresnel coefficient for P-polarised wave which is parallel to the plane of propogation/incidence. This function gives you back the reflectance. The transmittance is T=1-R */ + double Refl_P(double thetai); + + /* The temperature and attenuation model has been taken from AraSim which also took it from here http://icecube.wisc.edu/~araproject/radio/ . This is basically Matt Newcomb's icecube directory which has alot of information, plots and codes about South Pole Ice activities. Please read it if you find it interesting. */ + + /* Temperature model:The model takes in value of depth z in m and returns the value of temperature in Celsius.*/ + double GetIceTemperature(double z); + + /* Ice Attenuation Length model: Takes in value of frequency in Ghz and depth z and returns you the value of attenuation length in m */ + double GetIceAttenuationLength(double z, double frequency); + + /* Setup the integrand to calculate the attenuation */ + double AttenuationIntegrand (double x, void * params); + +/* Integrate over the integrand to calculate the attenuation */ + double IntegrateOverLAttn (double A0, double Frequency, double z0, double z1, double Lvalue); + +/* Calculate the total attenuation for each type of ray */ + double GetTotalAttenuationDirect (double A0, double frequency, double z0, double z1, double Lvalue); + + double GetTotalAttenuationReflected (double A0, double frequency, double z0, double z1, double Lvalue); + + double GetTotalAttenuationRefracted (double A0, double frequency, double z0, double z1, double zmax, double Lvalue); + + /* Use GSL minimiser which relies on calculating function deriavtives. This function uses GSL's Newton's algorithm to find root for a given function. */ + double FindFunctionRootFDF(gsl_function_fdf FDF,double x_lo, double x_hi); + + /* Use GSL minimiser which uses GSL's false position algorithm to find root for a given function. */ + double FindFunctionRoot(gsl_function F,double x_lo, double x_hi); + + /* Use GSL minimiser which uses GSL's false position algorithm to find root for a given function. */ + double FindFunctionRootZmax(gsl_function F,double x_lo, double x_hi); + + /* Define the function that will be minimized to get the value of the depth of the turning point for a given refracted ray. This function basically requires the value of the L parameter to find the depth of the turning point. This comes from the part of the fDnfR function where sqrt( n(z) - L ). This imposes the constraint then n(z)=L at the turning point of the ray from which we can find zmax. */ + struct Minnz_params { double a,l; }; + double GetMinnz(double x,void *params); + + /* Get the value of the depth of the turning point for the refracted ray */ + double GetZmax(double A, double L); + + /* Analytical solution describing ray paths in ice as function of depth */ + struct fDnfR_params { double a, b, c, l; }; + double fDnfR(double x,void *params); + + /* Analytical solution describing the ray path in ice as a function of the L parameter */ + struct fDnfR_L_params { double a, b, c, z; }; + double fDnfR_L(double x,void *params); + + /* The function used to calculate ray propogation time in ice */ + struct ftimeD_params { double a, b, c, speedc,l; }; + double ftimeD(double x,void *params); + + /* The function used to calculate ray geometric path in ice */ + double fpathD(double x,void *params); + + /* The set of functions starting with the name "fDa" are used in the minimisation procedure to find the launch angle (or the L parameter) for the direct ray */ + struct fDanfRa_params { double a, z0, x1, z1; }; + double fDa(double x,void *params); + + double fDa_df(double x,void *params); + + void fDa_fdf (double x, void *params,double *y, double *dy); + + /* The set of functions starting with the name "fRa" are used in the minimisation procedure to find the launch angle (or the L parameter) for the reflected ray */ + double fRa(double x,void *params); + + double fRa_df(double x,void *params); + + void fRa_fdf (double x, void *params,double *y, double *dy); + + /* This function is minimised to find the launch angle (or the L parameter) for the refracted ray */ + double fRaa(double x,void *params); + + double fRaa_df(double x,void *params); + + void fRaa_fdf (double x, void *params,double *y, double *dy); + + /* This functions works for the Direct ray and gives you back the launch angle, receive angle and propagation time of the ray together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. */ + double* GetDirectRayPar(double z0, double x1, double z1); + + /* This functions works for the Reflected ray and gives you back the launch angle, receive angle and propagation times (of the whole ray and the two direct rays that make it up) together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. */ + double *GetReflectedRayPar(double z0, double x1 ,double z1); + + /* This functions works for the Refracted ray and gives you back the launch angle, receive angle and propagation times (of the whole ray and the two direct rays that make it up) together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. It requires the launch angle of the reflected ray as an input. */ + double *GetRefractedRayPar(double z0, double x1 ,double z1, double LangR, double RangR, double checkzeroD, double checkzeroR); + + + /* This function returns the x and z values for the full Direct ray path and prints out the ray path in a text file */ + void GetFullDirectRayPath(double z0, double x1, double z1,double lvalueD, vector &x, vector &z); + + /* This function returns the x and z values for the full Reflected ray path and prints out the ray path in a text file */ + void GetFullReflectedRayPath(double z0, double x1, double z1,double lvalueR, vector &x, vector &z); + + /* This function returns the x and z values for the full Refracted ray path and prints out the ray path in a text file */ + void GetFullRefractedRayPath(double z0, double x1, double z1, double zmax, double lvalueRa, vector &x, vector &z,int raynumber); + + /* function for plotting and storing all the rays */ + void PlotAndStoreRays(double x0,double z0, double z1, double x1, double zmax[2], double lvalues[4], double checkzeroes[4]); + + double *IceRayTracing(double x0, double z0, double x1, double z1); + + /* Analytical solution describing ray paths in ice as function of depth for constant refractive index*/ + double fDnfR_Cnz(double x,void *params); + + /* Analytical solution describing the ray path in ice as a function of the L parameter for constant refractive index*/ + double fDnfR_L_Cnz(double x,void *params); + + /* This function is minimised to find the launch angle (or the L parameter) for the reflected ray for constant refractive index*/ + double fRa_Cnz(double x,void *params); + + double fRa_Cnz_df(double x,void *params); + + void fRa_Cnz_fdf (double x, void *params,double *y, double *dy); + + /* This functions works for the Direct ray and gives you back the launch angle, receive angle and propagation time of the ray together with values of the L parameter. This for constant refractive index*/ + double* GetDirectRayPar_Cnz(double z0, double x1, double z1, double A_ice_Cnz); + + /* This functions works for the Reflected ray and gives you back the launch angle, receive angle and propagation times (of the whole ray and the two direct rays that make it up) together with values of the L parameter. This is for constant refractive index*/ + double *GetReflectedRayPar_Cnz(double z0, double x1 , double z1, double A_ice_Cnz); + + /* This function returns the x and z values for the full Direct ray path and prints out the ray path in a text file. This is for a constant refractive index. */ + void GetFullDirectRayPath_Cnz(double z0, double x1, double z1, double lvalueD, double A_ice_Cnz,vector &x, vector &z); + + /* This function returns the x and z values for the full Reflected ray path and prints out the ray path in a text file. This is for a constant refractive index. */ + void GetFullReflectedRayPath_Cnz(double z0, double x1, double z1, double lvalueR, double A_ice_Cnz,vector &x, vector &z); + + /* function for plotting and storing all the rays. This is for constant refractive index. */ + void PlotAndStoreRays_Cnz(double x0,double z0, double z1, double x1, double lvalues[2], double A_ice_Cnz); + + /* This is the main raytracing function. x0 always has to be zero. z0 is the Tx depth in m and z1 is the depth of the Rx in m. Both depths are negative. x1 is the distance between them. This functions works for a constant refractive index */ + double *IceRayTracing_Cnz(double x0, double z0, double x1, double z1, double A_ice_Cnz); + + /* The set of functions starting with the name "fDa" are used in the minimisation procedure to find the launch angle (or the L parameter) for the direct ray */ + double fDa_Air(double x,void *params); + + /* This functions works for the Direct ray and gives you back the launch angle, receive angle and propagation time of the ray together with values of the L parameter and checkzero variable. checkzero variable checks how close the minimiser came to 0. 0 is perfect and less than 0.5 is pretty good. more than that should not be acceptable. */ + double* GetDirectRayPar_Air(double z0, double x1, double z1); + + double *DirectRayTracer(double xT, double yT, double zT, double xR, double yR, double zR); + + /* Function that makes interpolation tables for raytracing */ + void MakeTable(double ShowerHitDistance, double zT, int AntNum); + + /* Function that calculates the interpolated value for raytracing. The rt parameter: 0 is for D ray optical time, 1 is for D ray geometric path length, 2 is for D launch angle, 3 is for D recieve angle, 4 is D for ray attenuation, 5 is for R ray optical time, 6 is for R ray geometric path length, 7 is for R launch angle, 8 is for R recieve angle, 9 is R for ray attenuation */ + double GetInterpolatedValue(double xR, double zR, int rtParameter,int AntNum); + + //void GetRayTracingSolutions(double RxDepth, double Distance, double TxDepth, double TimeRay[2], double PathRay[2], double LaunchAngle[2], double RecieveAngle[2], int IgnoreCh[2], double IncidenceAngleInIce[2],vector xRay[2], vector zRay[2]); + void GetRayTracingSolutions(double RxDepth, double Distance, double TxDepth, double TimeRay[2], double PathRay[2], double LaunchAngle[2], double RecieveAngle[2], int IgnoreCh[2], double IncidenceAngleInIce[2], double A0, double frequency, double AttRay[2],vector xRay[2],vector zRay[2]); +} +#endif \ No newline at end of file diff --git a/Makefile b/Makefile index ce4ab3c..7431780 100644 --- a/Makefile +++ b/Makefile @@ -42,14 +42,14 @@ ARA_ROOT_FLAGS = # added for Fortran to C++ -LIBS = $(ROOTLIBS) -lMinuit $(SYSLIBS) -lRootFftwWrapper -lfftw3 +LIBS = $(ROOTLIBS) -lMinuit $(SYSLIBS) -lRootFftwWrapper -lfftw3 -lgsl -lgslcblas GLIBS = $(ROOTGLIBS) $(SYSLIBS) ROOT_LIBRARY = libAra.${DLLSUF} -OBJS = Vector.o EarthModel.o IceModel.o Birefringence.o Trigger.o Ray.o Tools.o Efficiencies.o Event.o Detector.o Position.o Spectra.o RayTrace.o RayTrace_IceModels.o signal.o secondaries.o Settings.o Primaries.o counting.o RaySolver.o Report.o eventSimDict.o AraSim.o -CCFILE = Vector.cc EarthModel.cc IceModel.cc Birefringence.cc Trigger.cc Ray.cc Tools.cc Efficiencies.cc Event.cc Detector.cc Spectra.cc Position.cc RayTrace.cc signal.cc secondaries.cc RayTrace_IceModels.cc Settings.cc Primaries.cc counting.cc RaySolver.cc Report.cc AraSim.cc -CLASS_HEADERS = Trigger.h Detector.h Settings.h Spectra.h IceModel.h Birefringence.h Primaries.h Report.h Event.h #need to add headers which added to Tree Branch +OBJS = Vector.o EarthModel.o IceModel.o Birefringence.o Trigger.o Ray.o Tools.o Efficiencies.o Event.o Detector.o Position.o Spectra.o RayTrace.o RayTrace_IceModels.o signal.o secondaries.o Settings.o Primaries.o counting.o RaySolver.o Report.o eventSimDict.o AraSim.o IceRayTracing.o +CCFILE = Vector.cc EarthModel.cc IceModel.cc Birefringence.cc Trigger.cc Ray.cc Tools.cc Efficiencies.cc Event.cc Detector.cc Spectra.cc Position.cc RayTrace.cc signal.cc secondaries.cc RayTrace_IceModels.cc Settings.cc Primaries.cc counting.cc RaySolver.cc Report.cc AraSim.cc IceRayTracing.cc +CLASS_HEADERS = Trigger.h Detector.h Settings.h Spectra.h IceModel.h Birefringence.h Primaries.h Report.h Event.h IceRayTracing.hh #need to add headers which added to Tree Branch PROGRAMS = AraSim @@ -99,7 +99,7 @@ endif eventSimDict.C: $(CLASS_HEADERS) @echo "Generating dictionary ..." @ rm -f *Dict* - rootcint -f $@ -c $(DICT_FLAGS) -I./ $(INC_ARA_UTIL) $(CLASS_HEADERS) ${ARA_ROOT_HEADERS} LinkDef.h + rootcint -f $@ -c $(DICT_FLAGS) -I./ $(INC_ARA_UTIL) $(SYSINCLUDES) $(CLASS_HEADERS) ${ARA_ROOT_HEADERS} LinkDef.h clean: @rm -f *Dict* diff --git a/RaySolver.cc b/RaySolver.cc index b734c0c..61dcba4 100644 --- a/RaySolver.cc +++ b/RaySolver.cc @@ -6,6 +6,7 @@ #include #include #include "Settings.h" +#include "IceRayTracing.hh" RaySolver::RaySolver() { @@ -1200,317 +1201,392 @@ void RaySolver::Solve_Ray (Position &source, Position &target, IceModel *antarct // test print out location //std::cout<<"--src_x="< 0 && sol_error > 0 ) { // this case do findPahts again with src, trg exchanged - paths_exc=tf.findPaths(trg,src,frequency/1.0e3,polarization, sol_cnt_exc, sol_error_exc, refl,requiredAccuracy); + double Distance= sqrt( (source_tmp.GetX()-target_tmp.GetX())*(source_tmp.GetX()-target_tmp.GetX()) + (source_tmp.GetY()-target_tmp.GetY())*(source_tmp.GetY()-target_tmp.GetY()) ) ; - std::cout<<"sol cnt exc : "< sol_error_exc ) { // exchanged one is better (less sol_error) - paths = paths_exc; - SrcTrgExc = 1; - } - else if (sol_error == sol_error_exc) { // sol_error are same for both paths. Then use the better miss dist. + IceRayTracing::GetRayTracingSolutions(source_tmp.GetZ(), Distance, target_tmp.GetZ(), TimeRay, PathRay, LaunchAngle, RecieveAngle, IgnoreCh, IncidenceAngleInIce, 1, 0.2,AttRay,xRay,zRay); + swap(LaunchAngle[0],RecieveAngle[0]); - std::vector::const_iterator it_tmp=paths.begin(); - std::vector::const_iterator it_tmp_exc=paths_exc.begin(); - if ( it_tmp->miss > it_tmp_exc->miss ) { // if fist one miss dist is worse - paths = paths_exc; - SrcTrgExc = 1; - } - } + TimeRay[0]=TimeRay[0];//*pow(10,-9); + TimeRay[1]=TimeRay[1];//*pow(10,-9); + + LaunchAngle[0]=(180-LaunchAngle[0])*(PI/180); + RecieveAngle[0]=(180-RecieveAngle[0])*(PI/180); + + LaunchAngle[1]=(LaunchAngle[1])*(PI/180); + RecieveAngle[1]=(RecieveAngle[1])*(PI/180); + + if(IncidenceAngleInIce[0]!=100){ + IncidenceAngleInIce[0]=IncidenceAngleInIce[0]*(PI/180); } - - - if(showLabels){ - if(!paths.empty()){ - //-------------------------------------------------- - // std::cout << std::fixed - // << "path length (m) " - // << "path time (ns) " - // << "launch angle " - // << "recipt angle " - // << "reflect angle " - // << "miss dist. " - // << "attenuation " - // << "amplitude" - // << std::endl; - //-------------------------------------------------- - solution_toggle = 1; - } - else { - //std::cout << "No solutions" << std::endl; - solution_toggle = 0; + if(IncidenceAngleInIce[1]!=100){ + IncidenceAngleInIce[1]=IncidenceAngleInIce[1]*(PI/180); + } + + solution_toggle = 0; + for(int iray=0;iray<2;iray++){ + if(IgnoreCh[iray]==1){ + outputs.resize(5); + + outputs[0].push_back(PathRay[iray]); + outputs[1].push_back(LaunchAngle[iray]); + outputs[2].push_back(RecieveAngle[iray]); + outputs[3].push_back(IncidenceAngleInIce[iray]); + outputs[4].push_back(TimeRay[iray]); // time in s (not ns) + + std::string pathfilename; + if ( sol_no == 0 ) { + pathfilename = "./pathfile_0.txt"; + } + else if ( sol_no == 1 ) { + pathfilename = "./pathfile_1.txt"; } - } - if(!dumpPaths){ - for(std::vector::const_iterator it=paths.begin(); it!=paths.end(); ++it){ - if ( it->sol_error == 0 ) { - //double signal = tf.signalStrength(*it,src,trg,refl); - //-------------------------------------------------- - // std::cout << std::left << std::fixed - // << std::setprecision(2) << std::setw(15) << it->pathLen << ' ' - // << std::setprecision(2) << std::setw(14) << 1e9*it->pathTime << ' ' - // << std::setprecision(4) << std::setw(12) << it->launchAngle << ' ' - // << std::setprecision(4) << std::setw(12) << it->receiptAngle << ' ' - // << std::setprecision(3) << std::setw(13) << it->reflectionAngle << ' ' - // << std::setprecision(2) << std::setw(10) << it->miss << ' ' - // << std::scientific << std::setprecision(4) << std::setw(11) << it->attenuation << ' ' - // //amplitude calculation, ignoring frequency response at both ends, angular response of receiver - // << std::setw(10) << (it->attenuation*signal) - // << std::endl; - //-------------------------------------------------- - if ( SrcTrgExc == 0 ) { // src, trg as original - outputs.resize(5); + solution_toggle = 1; + RayStep.resize(sol_no+1); - outputs[0].push_back(it->pathLen); - outputs[1].push_back(it->launchAngle); - outputs[2].push_back(it->receiptAngle); - outputs[3].push_back(it->reflectionAngle); - outputs[4].push_back( it->pathTime ); // time in s (not ns) - //std::cout<<"outputs[0]["<pathLen<<"\n"; + RayStep[sol_no].resize(2); // x and z for the direct + for(int step=0; steppathLen >= 10. ) // if more than 10m difference - { - //cout<<"source, target physical distance: "<pathLen<<", orgsol"<pathLen<<", orgsol"<pathLen<<", orgsol"<ANALYTIC_RAYTRACE_MODE==0){////do AraSim Numerical Raytracing + + RayTrace::TraceFinder tf(refractionModel,attenuationModel); + std::vector paths; + std::vector paths_exc; - std::string pathfilename; - if ( sol_no == 0 ) { - pathfilename = "./pathfile_0.txt"; - } - else if ( sol_no == 1 ) { - pathfilename = "./pathfile_1.txt"; - } - else { - pathfilename = "./pathfile.txt"; - } + paths=tf.findPaths(src,trg,frequency/1.0e3,polarization, sol_cnt, sol_error, refl,requiredAccuracy); - //testvector.resize(sol_no+1); - RayStep.resize(sol_no+1); + + //std::vector < std::vector < std::vector > > testvector; - // construct - pathStore_vector pathsave_test; - //pathStore_vector_2 pathsave_test (testvector); - //pathStore_test pathsave_test ("./pathfile.txt"); - //pathStore_test pathsave_test (pathfilename); - - //tf.doTrace(src_z, it->launchAngle, RayTrace::rayTargetRecord(trg_z,sqrt((trg_x-src_x)*(trg_x-src_x)+(trg_y-src_y)*(trg_y-src_y))), refl, 0.0, 0.0, sol_error, &pathsave_test); - tf.doTrace(src.GetZ(), it->launchAngle, RayTrace::rayTargetRecord(trg.GetZ(),sqrt((trg.GetX()-src.GetX())*(trg.GetX()-src.GetX())+(trg.GetY()-src.GetY())*(trg.GetY()-src.GetY()))), refl, 0.0, 0.0, sol_error, &pathsave_test); + //std::cout<<"sol cnt : "< 0 && sol_error > 0 ) { // this case do findPahts again with src, trg exchanged + paths_exc=tf.findPaths(trg,src,frequency/1.0e3,polarization, sol_cnt_exc, sol_error_exc, refl,requiredAccuracy); - //pathsave_test.CopyVector( testvector, sol_no ); - pathsave_test.CopyVector( RayStep, sol_no ); - pathsave_test.DelVector(); + std::cout<<"sol cnt exc : "< sol_error_exc ) { // exchanged one is better (less sol_error) + paths = paths_exc; + SrcTrgExc = 1; + } + else if (sol_error == sol_error_exc) { // sol_error are same for both paths. Then use the better miss dist. - //cout<<"\nRayStep : "<<(int)testvector[0].size()< 0 ) { - //dx = fabs(testvector[sol_no][0][step-1] - testvector[sol_no][0][step]); - //dz = fabs(testvector[sol_no][1][step-1] - testvector[sol_no][1][step]); - dx = fabs(RayStep[sol_no][0][step-1] - RayStep[sol_no][0][step]); - dz = fabs(RayStep[sol_no][1][step-1] - RayStep[sol_no][1][step]); - totalpath += sqrt( (dx*dx) + (dz*dz) ); + std::vector::const_iterator it_tmp=paths.begin(); + std::vector::const_iterator it_tmp_exc=paths_exc.begin(); + if ( it_tmp->miss > it_tmp_exc->miss ) { // if fist one miss dist is worse + paths = paths_exc; + SrcTrgExc = 1; + } + } + } + + + if(showLabels){ + if(!paths.empty()){ + //-------------------------------------------------- + // std::cout << std::fixed + // << "path length (m) " + // << "path time (ns) " + // << "launch angle " + // << "recipt angle " + // << "reflect angle " + // << "miss dist. " + // << "attenuation " + // << "amplitude" + // << std::endl; + //-------------------------------------------------- + solution_toggle = 1; + } + else { + //std::cout << "No solutions" << std::endl; + solution_toggle = 0; + } + } + if(!dumpPaths){ + for(std::vector::const_iterator it=paths.begin(); it!=paths.end(); ++it){ + if ( it->sol_error == 0 ) { + //double signal = tf.signalStrength(*it,src,trg,refl); + //-------------------------------------------------- + // std::cout << std::left << std::fixed + // << std::setprecision(2) << std::setw(15) << it->pathLen << ' ' + // << std::setprecision(2) << std::setw(14) << 1e9*it->pathTime << ' ' + // << std::setprecision(4) << std::setw(12) << it->launchAngle << ' ' + // << std::setprecision(4) << std::setw(12) << it->receiptAngle << ' ' + // << std::setprecision(3) << std::setw(13) << it->reflectionAngle << ' ' + // << std::setprecision(2) << std::setw(10) << it->miss << ' ' + // << std::scientific << std::setprecision(4) << std::setw(11) << it->attenuation << ' ' + // //amplitude calculation, ignoring frequency response at both ends, angular response of receiver + // << std::setw(10) << (it->attenuation*signal) + // << std::endl; + //-------------------------------------------------- + + if ( SrcTrgExc == 0 ) { // src, trg as original + outputs.resize(5); + + outputs[0].push_back(it->pathLen); + outputs[1].push_back(it->launchAngle); + outputs[2].push_back(it->receiptAngle); + outputs[3].push_back(it->reflectionAngle); + outputs[4].push_back( it->pathTime ); // time in s (not ns) + //std::cout<<"outputs[0]["<pathLen<<"\n"; + + + // test if raysol dist is shorter than physical distance + if ( distance_flat - it->pathLen >= 10. ) // if more than 10m difference + { + //cout<<"source, target physical distance: "<pathLen<<", orgsol"<pathLen<<", orgsol"<pathLen<<", orgsol"<pathLen<<", pathSum : "< pathsave_test; + //pathStore_vector_2 pathsave_test (testvector); + //pathStore_test pathsave_test ("./pathfile.txt"); + //pathStore_test pathsave_test (pathfilename); + + //tf.doTrace(src_z, it->launchAngle, RayTrace::rayTargetRecord(trg_z,sqrt((trg_x-src_x)*(trg_x-src_x)+(trg_y-src_y)*(trg_y-src_y))), refl, 0.0, 0.0, sol_error, &pathsave_test); + tf.doTrace(src.GetZ(), it->launchAngle, RayTrace::rayTargetRecord(trg.GetZ(),sqrt((trg.GetX()-src.GetX())*(trg.GetX()-src.GetX())+(trg.GetY()-src.GetY())*(trg.GetY()-src.GetY()))), refl, 0.0, 0.0, sol_error, &pathsave_test); + + //pathsave_test.CopyVector( testvector, sol_no ); + pathsave_test.CopyVector( RayStep, sol_no ); + pathsave_test.DelVector(); + + /* + double totalpath = 0.; + double dx, dz; + + //cout<<"\nRayStep : "<<(int)testvector[0].size()< 0 ) { + //dx = fabs(testvector[sol_no][0][step-1] - testvector[sol_no][0][step]); + //dz = fabs(testvector[sol_no][1][step-1] - testvector[sol_no][1][step]); + dx = fabs(RayStep[sol_no][0][step-1] - RayStep[sol_no][0][step]); + dz = fabs(RayStep[sol_no][1][step-1] - RayStep[sol_no][1][step]); + totalpath += sqrt( (dx*dx) + (dz*dz) ); + } + } + */ - //testvector.clear(); + //cout<<"pathLen : "<pathLen<<", pathSum : "<pathLen); - outputs[1].push_back(PI - it->receiptAngle); - outputs[2].push_back(PI - it->launchAngle); - outputs[3].push_back(it->reflectionAngle); - outputs[4].push_back( it->pathTime ); // time in s (not ns) - //std::cout<<"outputs[0]["<pathLen<<"\n"; + //testvector.clear(); - // test if raysol dist is shorter than physical distance - if ( distance_flat - it->pathLen >= 10. ) // if more than 10m difference - { - //cout<<"source, target physical distance: "<pathLen<<", excsol"<pathLen<<", orgsol"<pathLen<<", orgsol"<pathLen); + outputs[1].push_back(PI - it->receiptAngle); + outputs[2].push_back(PI - it->launchAngle); + outputs[3].push_back(it->reflectionAngle); + outputs[4].push_back( it->pathTime ); // time in s (not ns) + //std::cout<<"outputs[0]["<pathLen<<"\n"; + + + // test if raysol dist is shorter than physical distance + if ( distance_flat - it->pathLen >= 10. ) // if more than 10m difference + { + //cout<<"source, target physical distance: "<pathLen<<", excsol"<pathLen<<", orgsol"<pathLen<<", orgsol"< pathsave_test; + //pathStore_vector_2 pathsave_test (testvector); + //pathStore_test pathsave_test ("./pathfile.txt"); + //pathStore_test pathsave_test (pathfilename); + + //tf.doTrace(src_z, it->launchAngle, RayTrace::rayTargetRecord(trg_z,sqrt((trg_x-src_x)*(trg_x-src_x)+(trg_y-src_y)*(trg_y-src_y))), refl, 0.0, 0.0, sol_error, &pathsave_test); + tf.doTrace(src.GetZ(), PI - it->receiptAngle, RayTrace::rayTargetRecord(trg.GetZ(),sqrt((trg.GetX()-src.GetX())*(trg.GetX()-src.GetX())+(trg.GetY()-src.GetY())*(trg.GetY()-src.GetY()))), refl, 0.0, 0.0, sol_error, &pathsave_test); + + //pathsave_test.CopyVector( testvector, sol_no ); + pathsave_test.CopyVector( RayStep, sol_no ); + pathsave_test.DelVector(); + + /* + double totalpath = 0.; + double dx, dz; + + //cout<<"\nRayStep : "<<(int)testvector[0].size()< 0 ) { + //dx = fabs(testvector[sol_no][0][step-1] - testvector[sol_no][0][step]); + //dz = fabs(testvector[sol_no][1][step-1] - testvector[sol_no][1][step]); + dx = fabs(RayStep[sol_no][0][step-1] - RayStep[sol_no][0][step]); + dz = fabs(RayStep[sol_no][1][step-1] - RayStep[sol_no][1][step]); + totalpath += sqrt( (dx*dx) + (dz*dz) ); + } + } + */ - // construct - pathStore_vector pathsave_test; - //pathStore_vector_2 pathsave_test (testvector); - //pathStore_test pathsave_test ("./pathfile.txt"); - //pathStore_test pathsave_test (pathfilename); - - //tf.doTrace(src_z, it->launchAngle, RayTrace::rayTargetRecord(trg_z,sqrt((trg_x-src_x)*(trg_x-src_x)+(trg_y-src_y)*(trg_y-src_y))), refl, 0.0, 0.0, sol_error, &pathsave_test); - tf.doTrace(src.GetZ(), PI - it->receiptAngle, RayTrace::rayTargetRecord(trg.GetZ(),sqrt((trg.GetX()-src.GetX())*(trg.GetX()-src.GetX())+(trg.GetY()-src.GetY())*(trg.GetY()-src.GetY()))), refl, 0.0, 0.0, sol_error, &pathsave_test); + //cout<<"pathLen : "<pathLen<<", pathSum : "< 0 ) { - //dx = fabs(testvector[sol_no][0][step-1] - testvector[sol_no][0][step]); - //dz = fabs(testvector[sol_no][1][step-1] - testvector[sol_no][1][step]); - dx = fabs(RayStep[sol_no][0][step-1] - RayStep[sol_no][0][step]); - dz = fabs(RayStep[sol_no][1][step-1] - RayStep[sol_no][1][step]); - totalpath += sqrt( (dx*dx) + (dz*dz) ); - } + sol_no++; } - */ - //cout<<"pathLen : "<pathLen<<", pathSum : "<0){ + int num_solutions = outputs[0].size(); + if(num_solutions>1){ + // we only need to do something in the case of two solutions + double path_len_first_sol = outputs[0][0]; + double path_len_second_sol = outputs[0][1]; + if(path_len_second_sol < path_len_first_sol){ + printf("The second path length (%.2f) is shorter than the first (%.2f)\n",path_len_second_sol, path_len_first_sol); + printf("Flipping the order in the output vectors..."); + + // if the shorter path length (which must be the direct solution) + // has ended up second, then we need to make a swap + // so that it's first in the output + + // for the outputs vector, we flip the first and second columns + // that is, the direct and refracted/reflected solution + // first, make a copy of the original + vector< vector > tmp_outputs(outputs); + // then make the flip + for(int item=0; item > > tmp_RayStep(RayStep); + // now, totally empty the older vector (the one we need to return to the user at the end) + RayStep.clear(); + RayStep.resize(2); // the RayStep needs a direct and refr/refl solution + // first, the direct solution (which is currently in tmp_RayStep[1]) + RayStep[0].resize(2); // x and z for the direct + for(int step=0; step print; + for(std::vector::const_iterator it=paths.begin(); it!=paths.end(); ++it){ + tf.doTrace(src_z, it->launchAngle, RayTrace::rayTargetRecord(trg_z,sqrt((trg_x-src_x)*(trg_x-src_x)+(trg_y-src_y)*(trg_y-src_y))), refl, 0.0, 0.0, sol_error, &print); + //std::cout << "\n\n"; + } + } + //-------------------------------------------------- + // return(0); + //-------------------------------------------------- - } - - // reorder the solutions (BAC Jan 2021) - // it is useful at this point to re-order the solutions so that the direct solution - // (if it exists) is reported first - // this copying method is rather inefficient in terms of memory allocation - // but also, it's not obvious that it ever gets called - // (Brian couldn't find a case in testing on a few hundred events) - // so probably it doesn't matter - if (outputs.size()>0){ - int num_solutions = outputs[0].size(); - if(num_solutions>1){ - // we only need to do something in the case of two solutions - double path_len_first_sol = outputs[0][0]; - double path_len_second_sol = outputs[0][1]; - if(path_len_second_sol < path_len_first_sol){ - printf("The second path length (%.2f) is shorter than the first (%.2f)\n",path_len_second_sol, path_len_first_sol); - printf("Flipping the order in the output vectors..."); - - // if the shorter path length (which must be the direct solution) - // has ended up second, then we need to make a swap - // so that it's first in the output - - // for the outputs vector, we flip the first and second columns - // that is, the direct and refracted/reflected solution - // first, make a copy of the original - vector< vector > tmp_outputs(outputs); - // then make the flip - for(int item=0; item > > tmp_RayStep(RayStep); - - // now, totally empty the older vector (the one we need to return to the user at the end) - RayStep.clear(); - RayStep.resize(2); // the RayStep needs a direct and refr/refl solution - - // first, the direct solution (which is currently in tmp_RayStep[1]) - RayStep[0].resize(2); // x and z for the direct - for(int step=0; step print; - for(std::vector::const_iterator it=paths.begin(); it!=paths.end(); ++it){ - tf.doTrace(src_z, it->launchAngle, RayTrace::rayTargetRecord(trg_z,sqrt((trg_x-src_x)*(trg_x-src_x)+(trg_y-src_y)*(trg_y-src_y))), refl, 0.0, 0.0, sol_error, &print); - //std::cout << "\n\n"; - } - } - //-------------------------------------------------- - // return(0); - //-------------------------------------------------- } diff --git a/Settings.cc b/Settings.cc index 60bc966..1b9ee5b 100644 --- a/Settings.cc +++ b/Settings.cc @@ -294,6 +294,8 @@ outputdir="outputs"; // directory where outputs go RAY_TRACE_ICE_MODEL_PARAMS=0; // Default: South Pole values fitted from RICE data + ANALYTIC_RAYTRACE_MODE=0; //default: 0 -- use numerical RayTracing for AraSim, 1 -- use analytical raytracing + WAVEFORM_LENGTH = 64/2*20; // Default: 64 digitization samples per block / 2 samples per waveform value * 20 blocks (value used for 2013-2016) WAVEFORM_CENTER = 0; // Default: 0, no offset in waveform centering @@ -318,6 +320,7 @@ outputdir="outputs"; // directory where outputs go APPLY_NOISE_FIGURE=0; // default: 0 - don't use new noise figure information CUSTOM_ELECTRONICS=0; //default: 0 -- don't use custom electronics, load regular "ARA_Electronics_TotalGain_TwoFilter.csv" + ELECTRONICS_ANTENNA_CONSISTENCY = 1; // default: 1 -- ensure antenna gain used to calculate electronics gain is consistent // with that used in this simulation @@ -694,6 +697,9 @@ void Settings::ReadFile(string setupfile) { else if (label == "RAY_TRACE_ICE_MODEL_PARAMS") { RAY_TRACE_ICE_MODEL_PARAMS = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); } + else if (label == "ANALYTIC_RAYTRACE_MODE"){ + ANALYTIC_RAYTRACE_MODE = atoi(line.substr(line.find_first_of("=") + 1).c_str()); + } else if (label == "WAVEFORM_LENGTH") { WAVEFORM_LENGTH = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); } diff --git a/Settings.h b/Settings.h index cf46fec..61bbe02 100644 --- a/Settings.h +++ b/Settings.h @@ -295,6 +295,8 @@ class Settings //30 : Mizuho (Ebimuna (1983)) //40: UNL Modified (PA model). Related slide: https://aradocs.wipac.wisc.edu/docs/0022/002222/001/inIceMC_Hughes_A5locations_10222020.pdf + int ANALYTIC_RAYTRACE_MODE; //default: 0 -- use numerical RayTracing for AraSim, 1 -- use analytical raytracing + int WAVEFORM_LENGTH; // the number of samples in the waveform length for V_mimic and UsefulAtriStationEvent, default: 64/2*20 = 640 int WAVEFORM_CENTER; // the relative location of the center of the write-out window with respect to the last triggered bin (which is laced at the center of the window by default), this effectively provides a global delay in the write-out window across all channels: positive values shift the write-out window to later times in the waveform, negative values shift the window to earlier times, default: 0