Computes the spin of the final BH that gets formed after merger. This is done usingn Eq 5-6 of arXiv:0710.3345
321058
def FinalSpin( Xi, eta ):
    """"""
    s4 = -0.129
    s5 = -0.384
    t0 = -2.686
    t2 = -3.454
    t3 = 2.353
    etaXi = eta * Xi
    eta2 = eta*eta
    finspin = (Xi + s4*Xi*etaXi + s5*etaXi*eta + t0*etaXi + 2.*(3.**0.5)*eta + t2*eta2 + t3*eta2*eta)
    if finspin > 1.0:
        raise ValueError("Value of final spin > 1.0. Aborting")
    else:
        return finspin
929826
def phi_a(mass1, mass2, spin1x, spin1y, spin2x, spin2y):
    
    phi1 = phi_from_spinx_spiny(primary_spin(mass1, mass2, spin1x, spin2x),
                                primary_spin(mass1, mass2, spin1y, spin2y))
    phi2 = phi_from_spinx_spiny(secondary_spin(mass1, mass2, spin1x, spin2x),
                                secondary_spin(mass1, mass2, spin1y, spin2y))
    return (phi1 - phi2) % (2 * numpy.pi)
953608
def secondary_xi(mass1, mass2, spin1x, spin1y, spin2x, spin2y):
    
    spinx = secondary_spin(mass1, mass2, spin1x, spin2x)
    spiny = secondary_spin(mass1, mass2, spin1y, spin2y)
    return xi2_from_mass1_mass2_spin2x_spin2y(mass1, mass2, spinx, spiny)
316704
def phi_from_spinx_spiny(spinx, spiny):
    
    phi = numpy.arctan2(spiny, spinx)
    return phi % (2 * numpy.pi)
878558
def _proc_no_rot_sym(self):
        
        self.sch_symbol = "C1"
        if self.is_valid_op(PointGroupAnalyzer.inversion_op):
            self.sch_symbol = "Ci"
            self.symmops.append(PointGroupAnalyzer.inversion_op)
        else:
            for v in self.principal_axes:
                mirror_type = self._find_mirror(v)
                if not mirror_type == "":
                    self.sch_symbol = "Cs"
                    break
473464
def canonical(self):
        
        Ac1 = np.hstack((self.deltah, np.zeros((self.nh, self.nz))))
        Ac2 = np.hstack((np.zeros((self.nz, self.nh)), self.a22))
        Ac = np.vstack((Ac1, Ac2))
        Bc = np.vstack((self.thetah, np.zeros((self.nz, self.nc))))
        Cc = np.vstack((np.zeros((self.nh, self.nw)), self.c2))
        Rc1 = np.hstack((self.llambda.T.dot(self.llambda), -
                         self.llambda.T.dot(self.ub)))
        Rc2 = np.hstack((-self.ub.T.dot(self.llambda), self.ub.T.dot(self.ub)))
        Rc = np.vstack((Rc1, Rc2))
        Qc = self.pih.T.dot(self.pih)
        Nc = np.hstack(
            (self.pih.T.dot(self.llambda), -self.pih.T.dot(self.ub)))

        lq_aux = LQ(Qc, Rc, Ac, Bc, N=Nc, beta=self.beta)

        P1, F1, d1 = lq_aux.stationary_values()

        self.F_b = F1[:, 0:self.nh]
        self.F_f = F1[:, self.nh:]

        self.pihat = np.linalg.cholesky(self.pih.T.dot(
            self.pih) + self.beta.dot(self.thetah.T).dot(P1[0:self.nh, 0:self.nh]).dot(self.thetah)).T
        self.llambdahat = self.pihat.dot(self.F_b)
        self.ubhat = - self.pihat.dot(self.F_f)

        return
979174
private void calculateRotationalAxis(Matrix rotation, Atom translation,
			double c) {
		// Calculate magnitude of rotationAxis components, but not signs
		double sum=0;
		double[] rotAx = new double[3];
		for(int i=0;i<3;i++) {
			rotAx[i] = Math.sqrt(rotation.get(i, i)-c);
			sum+=rotAx[i]*rotAx[i];
		}
		for(int i=0;i<3;i++) {
			rotAx[i] /= Math.sqrt(sum);
		}

		// Now determine signs
		double d0 = rotation.get(2,1)-rotation.get(1,2); //=2u[0]*sin(theta)
		double d1 = rotation.get(0,2)-rotation.get(2,0); //=2u[1]*sin(theta)
		double d2 = rotation.get(1,0)-rotation.get(0,1); //=2u[2]*sin(theta)

		double s12 = rotation.get(2,1)+rotation.get(1,2); //=2*u[1]*u[2]*(1-cos(theta))
		double s02 = rotation.get(0,2)+rotation.get(2,0); //=2*u[0]*u[2]*(1-cos(theta))
		double s01 = rotation.get(1,0)+rotation.get(0,1); //=2*u[0]*u[1]*(1-cos(theta))

		//Only use biggest d for the sign directly, for numerical stability around 180deg
		if( Math.abs(d0) < Math.abs(d1) ) { // not d0
			if( Math.abs(d1) < Math.abs(d2) ) { //d2
				if(d2>=0){ //u[2] positive
					if( s02 < 0 ) rotAx[0] = -rotAx[0];
					if( s12 < 0 ) rotAx[1] = -rotAx[1];
				} else { //u[2] negative
					rotAx[2] = -rotAx[2];
					if( s02 >= 0 ) rotAx[0] = -rotAx[0];
					if( s12 >= 0 ) rotAx[1] = -rotAx[1];
				}
			} else { //d1
				if(d1>=0) {//u[1] positive
					if( s01 < 0) rotAx[0] = -rotAx[0];
					if( s12 < 0) rotAx[2] = -rotAx[2];
				} else { //u[1] negative
					rotAx[1] = -rotAx[1];
					if( s01 >= 0) rotAx[0] = -rotAx[0];
					if( s12 >= 0) rotAx[2] = -rotAx[2];
				}
			}
		} else { // not d1
			if( Math.abs(d0) < Math.abs(d2) ) { //d2
				if(d2>=0){ //u[2] positive
					if( s02 < 0 ) rotAx[0] = -rotAx[0];
					if( s12 < 0 ) rotAx[1] = -rotAx[1];
				} else { //u[2] negative
					rotAx[2] = -rotAx[2];
					if( s02 >= 0 ) rotAx[0] = -rotAx[0];
					if( s12 >= 0 ) rotAx[1] = -rotAx[1];
				}
			} else { //d0
				if(d0>=0) { //u[0] positive
					if( s01 < 0 ) rotAx[1] = -rotAx[1];
					if( s02 < 0 ) rotAx[2] = -rotAx[2];
				} else { //u[0] negative
					rotAx[0] = -rotAx[0];
					if( s01 >= 0 ) rotAx[1] = -rotAx[1];
					if( s02 >= 0 ) rotAx[2] = -rotAx[2];
				}
			}
		}

		rotationAxis = new AtomImpl();
		rotationAxis.setCoords(rotAx);

		// Calculate screw = (rotationAxis dot translation)*u
		double dotProduct = Calc.scalarProduct(rotationAxis, translation);

		screwTranslation = Calc.scale(rotationAxis, dotProduct);
		otherTranslation = Calc.subtract(translation, screwTranslation);

		Atom hypot = Calc.vectorProduct(otherTranslation,rotationAxis);
		Calc.scaleEquals(hypot,.5/Math.tan(theta/2.0));

		// Calculate rotation axis position
		rotationPos = Calc.scaleAdd(.5,otherTranslation, hypot);
	}
395359
def coreBurkAlpha(self, R, Rs, rho0, r_core, ax_x, ax_y):
        
        x = R * Rs ** -1
        p = Rs * r_core ** -1

        gx = self._G(x, p)

        a = 2 * rho0 * Rs ** 2 * gx / x

        return a * ax_x / R, a * ax_y / R
1358660
def calculate_angular_momentum(self):
        
        clibrebound.reb_tools_angular_momentum.restype = reb_vec3d
        L = clibrebound.reb_tools_angular_momentum(byref(self))
        return [L.x, L.y, L.z]
333383
def spin2_polar(self):
        """"""
        return coordinates.cartesian_to_spherical_polar(
                                     self.spin2x, self.spin2y, self.spin2z)
873883
private static double combineRmsd(int b1, int b2, AFPChain afpChain,Atom[] ca1,Atom[] ca2)
	{
		int     i;
		int     afpn = 0;

		int[] afpChainList =afpChain.getAfpChainList();

		int[] block2Afp = afpChain.getBlock2Afp();
		int[] blockSize = afpChain.getBlockSize();


		int[]   list = new int[blockSize[b1]+blockSize[b2]];
		for(i = block2Afp[b1]; i < block2Afp[b1] + blockSize[b1]; i ++) {
			list[afpn ++] = afpChainList[i];
		}
		for(i = block2Afp[b2]; i < block2Afp[b2] + blockSize[b2]; i ++) {
			list[afpn ++] = afpChainList[i];
		}
		double  rmsd = AFPChainer.calAfpRmsd(afpn, list,0, afpChain,ca1,ca2);

		afpChain.setBlock2Afp(block2Afp);
		afpChain.setBlockSize(blockSize);
		afpChain.setAfpChainList(afpChainList);

		return rmsd;
	}
360338
def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1):
    
    implemented_phase_order = 7
    implemented_spin_order = 7
    if phase_order > implemented_phase_order:
        raise ValueError("pN coeffiecients of that order have not been implemented")
    elif phase_order == -1:
        phase_order = implemented_phase_order

    if spin_order > implemented_spin_order:
        raise ValueError("pN coeffiecients of that order have not been implemented")
    elif spin_order == -1:
        spin_order = implemented_spin_order

    qmdef1 = 1.0
    qmdef2 = 1.0

    M = m1 + m2
    dm = (m1-m2)/M
    m1M = m1 / M
    m2M = m2 / M

    s1z = s1z * m1M * m1M
    s2z = s2z * m2M * m2M

    _, eta = mass1_mass2_to_mchirp_eta(m1, m2)

    ecof = numpy.zeros(phase_order+1)
    # Orbital terms
    if phase_order >= 0:
        ecof[0] = 1.0
    if phase_order >= 1:
        ecof[1] = 0
    if phase_order >= 2:
        ecof[2] = -(1.0/12.0) * (9.0 + eta)
    if phase_order >= 3:
        ecof[3] = 0
    if phase_order >= 4:
        ecof[4] = (-81.0 + 57.0*eta - eta*eta) / 24.0
    if phase_order >= 5:
        ecof[5] = 0
    if phase_order >= 6:
        ecof[6] = - 675.0/64.0 + ( 34445.0/576.0    \
              - 205.0/96.0 * lal.PI * lal.PI ) * eta  \
              - (155.0/96.0) *eta * eta - 35.0/5184.0 * eta * eta
    # Spin terms

    ESO15s1 = 8.0/3.0 + 2.0*m2/m1
    ESO15s2 = 8.0/3.0 + 2.0*m1/m2

    ESS2 = 1.0 / eta
    EQM2s1 = qmdef1/2.0/m1M/m1M
    EQM2s1L = -qmdef1*3.0/2.0/m1M/m1M
    #EQM2s2 = qmdef2/2.0/m2M/m2M
    EQM2s2L = -qmdef2*3.0/2.0/m2M/m2M

    ESO25s1 = 11.0 - 61.0*eta/9.0 + (dm/m1M) * (-3.0 + 10.*eta/3.0)
    ESO25s2 = 11.0 - 61.0*eta/9.0 + (dm/m2M) * (3.0 - 10.*eta/3.0)

    ESO35s1 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 + (dm/m1M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0)
    ESO35s2 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 - (dm/m2M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0)

    if spin_order >=3:
        ecof[3] += ESO15s1 * s1z + ESO15s2 * s2z
    if spin_order >=4:
        ecof[4] += ESS2 * (s1z*s2z - 3.0*s1z*s2z)
        ecof[4] += EQM2s1*s1z*s1z + EQM2s1*s2z*s2z + EQM2s1L*s1z*s1z + EQM2s2L*s2z*s2z
    if spin_order >=5:
        ecof[5] = ESO25s1*s1z + ESO25s2*s2z
    if spin_order >=7:
        ecof[7] += ESO35s1*s1z + ESO35s2*s2z

    return ecof
1269780
def _alpha2rho0(self, theta_Rs, Rs):

        

        rho0 = theta_Rs / (4. * Rs ** 2 * (1. + np.log(1. / 2.)))
        return rho0
723898
private static boolean isLeftTurn(Point2D a, Point2D b, Point2D c) {
        return winding(a, b, c) > 0;
    }
965539
def spher_harms(l, m, inclination):
    

    # FIXME: we are using spin -2 weighted spherical harmonics for now,
    # when possible switch to spheroidal harmonics.
    Y_lm = lal.SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m).real
    Y_lminusm = lal.SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m).real
    Y_plus = Y_lm + (-1)**l * Y_lminusm
    Y_cross = Y_lm - (-1)**l * Y_lminusm

    return Y_plus, Y_cross
345021
def spintaylorf2(**kwds):
    
    #####Pull out the input arguments#####
    f_lower = double(kwds['f_lower'])
    delta_f = double(kwds['delta_f'])
    distance = double(kwds['distance'])
    mass1 = double(kwds['mass1'])
    mass2 = double(kwds['mass2'])
    spin1x = double(kwds['spin1x'])
    spin1y = double(kwds['spin1y'])
    spin1z = double(kwds['spin1z'])
    phi0 = double(kwds['coa_phase'])               #Orbital Phase at coalescence
    phase_order = int(kwds['phase_order'])
    amplitude_order = int(kwds['amplitude_order'])
    inclination = double(kwds['inclination'])
    lnhatx = sin(inclination)
    lnhaty = 0.
    lnhatz = cos(inclination)
    psi = 0.

    tC= -1.0 / delta_f
    M = mass1 + mass2
    eta = mass1 * mass2 / (M * M)
    m_sec = M * lal.MTSUN_SI
    piM = lal.PI * m_sec

    vISCO = 1. / sqrt(6.)
    fISCO = vISCO * vISCO * vISCO / piM
    f_max = ceilpow2(fISCO)
    n = int(f_max / delta_f + 1)
    kmax = int(fISCO / delta_f)
    kmin = int(numpy.ceil(f_lower / delta_f))
    kmax = kmax if (kmax<n) else n

    #####Calculate the Orientation#####
    v0 = pow(piM *  kmin * delta_f,1./3)
    chi = sqrt(spin1x**2+spin1y**2+spin1z**2)
    kappa = (lnhatx*spin1x+lnhaty*spin1y+lnhatz*spin1z)/chi if (chi > 0.)  else 1.
    Jx0 = mass1*mass2*lnhatx/v0 + mass1*mass1*spin1x
    Jy0 = mass1*mass2*lnhaty/v0 + mass1*mass1*spin1y
    Jz0 = mass1*mass2*lnhatz/v0 + mass1*mass1*spin1z
    thetaJ = acos(Jz0 / sqrt(Jx0**2+Jy0**2+Jz0**2))
    psiJ = atan2(Jy0, -Jx0) # FIXME: check that Jy0 and Jx0 are not both 0
    # Rotate Lnhat back to frame where J is along z, to figure out initial alpha
    rotLx = lnhatx*cos(thetaJ)*cos(psiJ) - lnhaty*cos(thetaJ)*sin(psiJ) + lnhatz*sin(thetaJ)
    rotLy = lnhatx*sin(psiJ) + lnhaty*cos(psiJ)
    alpha0 = atan2(rotLy, rotLx) # FIXME: check that rotLy and rotLx are not both 0
    psiJ_P =psiJ + psi
    psiJ_C =psiJ + psi + lal.PI/4.

    #####Calculate the Coefficients#####
    #quadparam = 1.
    gamma0 = mass1*chi/mass2
    #Calculate the spin corrections
    # FIXME should use pycbc's function, but sigma has different expression
    # in Andy's code, double check
    # pn_beta, pn_sigma, pn_gamma = pycbc.pnutils.mass1_mass2_spin1z_spin2z_to_beta_sigma_gamma(
    #                               mass1, mass2, chi*kappa, 0) # FIXME: spin2 is taken to be 0
    pn_beta = (113.*mass1/(12.*M) - 19.*eta/6.)*chi*kappa
    pn_sigma = (  (5.*(3.*kappa*kappa-1.)/2.) + (7. - kappa*kappa)/96.  ) * (mass1*mass1*chi*chi/M/M)
    pn_gamma = (5.*(146597. + 7056.*eta)*mass1/(2268.*M) - 10.*eta*(1276. + 153.*eta)/81.)*chi*kappa
    prec_fac0 = 5.*(4. + 3.*mass2/mass1)/64.
    dtdv2 = 743./336. + 11.*eta/4.
    dtdv3 = -4.*lal.PI + pn_beta
    dtdv4 = 3058673./1016064. + 5429.*eta/1008. + 617.*eta*eta/144. - pn_sigma
    dtdv5 = (-7729./672.+13.*eta/8.)*lal.PI + 9.*pn_gamma/40.

    #####Calculate the Initial Euler Angles alpha_ref, beta_ref=0 and zeta_ref#####
    gam = gamma0*v0
    sqrtfac = sqrt(1. + 2.*kappa*gam + gam*gam)
    logv0 = log(v0)
    logfac1 = log(1. + kappa*gam + sqrtfac)
    logfac2 = log(kappa + gam + sqrtfac)
    v02 = v0 * v0
    v03 = v0 * v02
    kappa2 = kappa * kappa
    kappa3 = kappa2 * kappa
    gamma02 = gamma0 * gamma0
    gamma03 = gamma02 *gamma0

    alpha_ref = prec_fac0*(  logfac2 *( dtdv2*gamma0 + dtdv3*kappa - dtdv5*kappa/(2.*gamma02) + dtdv4/(2.*gamma0) - dtdv4*kappa2/(2.*gamma0) + (dtdv5*kappa3)/(2.*gamma02) )  +  logfac1*( - dtdv2*gamma0*kappa - dtdv3 + kappa*gamma03/2. - gamma03*kappa3/2. ) + logv0 *( dtdv2*gamma0*kappa + dtdv3 - kappa*gamma03/2. + gamma03*kappa3/2. ) + sqrtfac *( dtdv3 + dtdv4*v0/2. + dtdv5/gamma02/3. + dtdv4*kappa/(2.*gamma0) + dtdv5*kappa*v0/(6.*gamma0) - dtdv5*kappa2/(2.*gamma02) - 1/(3.*v03) - gamma0*kappa/(6.*v02) - dtdv2/v0 - gamma02/(3.*v0) + gamma02*kappa2/(2.*v0) + dtdv5*v02/3. ))  - alpha0

    zeta_ref = prec_fac0*( dtdv3*gamma0*kappa*v0 + dtdv4*v0 + logfac2 *(-dtdv2*gamma0 - dtdv3*kappa + dtdv5*kappa/(2.*gamma02) - dtdv4/(2.*gamma0) + dtdv4*kappa2/(2.*gamma0) - dtdv5*kappa3/(2.*gamma02) ) + logv0 *( kappa*gamma03/2. - gamma03*kappa3/2. ) + logfac1 *( dtdv2*gamma0*kappa + dtdv3 - kappa*gamma03/2. + gamma03*kappa3/2. ) - 1/(3.*v03) - gamma0*kappa/(2.*v02) - dtdv2/v0 + dtdv4*gamma0*kappa*v02/2. + dtdv5*v02/2. + sqrtfac *( -dtdv3 - dtdv4*v0/2. - dtdv5/(3.*gamma02) - dtdv4*kappa/(2.*gamma0) - dtdv5*kappa*v0/(6.*gamma0) + dtdv5*kappa2/(2.*gamma02) + 1/(3.*v03) + gamma0*kappa/(6.*v02) + dtdv2/v0 + gamma02/(3.*v0) - gamma02*kappa2/(2.*v0) - dtdv5*v02/3. ) + dtdv5*gamma0*kappa*v03/3. )

    #####Calculate the Complex sideband factors, mm=2 is first entry#####
    RE_SBfac0= (1.+cos(thetaJ)**2)/2.
    RE_SBfac1= sin(2.*thetaJ)
    RE_SBfac2= 3.*sin(thetaJ)**2
    RE_SBfac3= -sin(2.*thetaJ)
    RE_SBfac4= (1.+cos(thetaJ)**2)/2.
    IM_SBfac0= -cos(thetaJ)
    IM_SBfac1= -2.*sin(thetaJ)
    IM_SBfac2= 0.
    IM_SBfac3= -2.*sin(thetaJ)
    IM_SBfac4= cos(thetaJ)

    #####Calculate the PN terms # FIXME replace with functions in lalsimulation #####
    theta = -11831./9240.
    lambdaa = -1987./3080.0
    pfaN = 3.0/(128.0 * eta)
    pfa2 = 5.0*(743.0/84 + 11.0 * eta)/9.0
    pfa3 = -16.0*lal.PI + 4.0*pn_beta
    pfa4 = 5.0*(3058.673/7.056 + 5429.0/7.0 * eta + 617.0 * eta*eta)/72.0 - \
            10.0*pn_sigma
    pfa5 = 5.0/9.0 * (7729.0/84.0 - 13.0 * eta) * lal.PI - pn_gamma
    pfl5 = 5.0/3.0 * (7729.0/84.0 - 13.0 * eta) * lal.PI - pn_gamma * 3
    pfa6 = (11583.231236531/4.694215680 - 640.0/3.0 * lal.PI * lal.PI- \
            6848.0/21.0*lal.GAMMA) + \
            eta * (-15335.597827/3.048192 + 2255./12. * lal.PI * \
            lal.PI - 1760./3.*theta +12320./9.*lambdaa) + \
            eta*eta * 76055.0/1728.0 - \
            eta*eta*eta*  127825.0/1296.0
    pfl6 = -6848.0/21.0
    pfa7 = lal.PI * 5.0/756.0 * ( 15419335.0/336.0 + 75703.0/2.0 * eta - \
            14809.0 * eta*eta)

    FTaN = 32.0 * eta*eta / 5.0
    FTa2 = -(12.47/3.36 + 3.5/1.2 * eta)
    FTa3 = 4.0 * lal.PI
    FTa4 = -(44.711/9.072 - 92.71/5.04 * eta - 6.5/1.8 * eta*eta)
    FTa5 = -(81.91/6.72 + 58.3/2.4 * eta) * lal.PI
    FTa6 = (664.3739519/6.9854400 + 16.0/3.0 * lal.PI*lal.PI -
            17.12/1.05 * lal.GAMMA +
         (4.1/4.8 * lal.PI*lal.PI - 134.543/7.776) * eta -
         94.403/3.024 * eta*eta - 7.75/3.24 * eta*eta*eta)
    FTl6 = -8.56/1.05
    FTa7 = -(162.85/5.04 - 214.745/1.728 * eta - 193.385/3.024 * eta*eta) \
            * lal.PI

    dETaN = 2 * -eta/2.0
    dETa1 = 2 * -(3.0/4.0 + 1.0/12.0 * eta)
    dETa2 = 3 * -(27.0/8.0 - 19.0/8.0 * eta + 1./24.0 * eta*eta)
    dETa3 = 4 * -(67.5/6.4 - (344.45/5.76 - 20.5/9.6 * lal.PI*lal.PI) *
                             eta + 15.5/9.6 * eta*eta + 3.5/518.4 * eta*eta*eta)

    amp0 = -4. * mass1 * mass2 / (1.0e+06 * distance * lal.PC_SI ) * \
                    lal.MRSUN_SI * lal.MTSUN_SI * sqrt(lal.PI/12.0)

    htildeP = FrequencySeries(zeros(n,dtype=complex128), delta_f=delta_f, copy=False)
    htildeC = FrequencySeries(zeros(n,dtype=complex128), delta_f=delta_f, copy=False)
    spintaylorf2_kernel(htildeP.data[kmin:kmax], htildeC.data[kmin:kmax],
                        kmin, phase_order, amplitude_order, delta_f, piM, pfaN,
                        pfa2, pfa3, pfa4, pfa5, pfl5,
                        pfa6, pfl6, pfa7, FTaN, FTa2,
                        FTa3, FTa4, FTa5, FTa6,
                        FTl6, FTa7, dETaN, dETa1, dETa2,  dETa3,
                        amp0, tC, phi0,
                        kappa, prec_fac0, alpha_ref, zeta_ref,
                        dtdv2, dtdv3, dtdv4, dtdv5,
                        RE_SBfac0, RE_SBfac1, RE_SBfac2, RE_SBfac3, RE_SBfac4,
                        IM_SBfac0, IM_SBfac1, IM_SBfac2, IM_SBfac3, IM_SBfac4,
                        psiJ_P, psiJ_C, gamma0)

    return htildeP, htildeC
772870
def _get_b(self):
        
        orig = self.origin
        dest = self.destination
        (ox, oy) = orig.pos
        (dx, dy) = dest.pos
        denominator = dx - ox
        x_numerator = (dy - oy) * ox
        y_numerator = denominator * oy
        return ((y_numerator - x_numerator), denominator)