MCMAC.py

'''
Version 0.0
- The initial version of the old TSM_v11 that is now being tracked in Git.
'''
from __future__ import division
import numpy
import scipy.integrate
import pickle
import profiles
import time
import sys
import cosmo

# Constants and conversions
G = 4.3*10**(-9) #newton's constant in units of Mpc*(km/s)**2 / M_sun
c = 3e5 #speed of light km/s
sinGyr = 31556926.*10**9 # s in a Giga-year
kginMsun = 1.98892*10**30 # kg in a solar mass
kminMpc = 3.08568025*10**19 # km in a Megaparsec
minMpc = 3.08568025*10**22 # m in a Megaparsec            
r2d = 180/numpy.pi # radians to degrees conversion factor

def randdraw(A,index=None):
    '''
    This function randomly draws a variable from either and input gaussian
    distribution or an array of random draws from a parent distribution function
    (e.g. bootstrap sampling).
    Input:
    A = [(float, float) or (1D array of floats)] If (float, float) is input then
        it is interpreted as (mu, sigma) of a Gaussian distribution. If (1D
        array of floats with length > 2) then it is interpreted as a random
        sample from a parent distribution.
    index = [int] rather than drawing a random value from the array the value
    Output:
    a = [float] a random draw from either the Gaussian or sample distribution.
    '''
    N_A = numpy.size(A)
    if N_A == 2: #then assumed (mu, sigma) format
        a = A[0]+A[1]*numpy.random.randn()
    elif N_A > 2: #Then assume distribution array format
        if index == None:
            index = numpy.random.randint(N_A)
        a = A[index]
    elif N_A <2:
        print 'MCMAC.randdraw: Error, parameter input array is not a valid size, exiting'
        sys.exit()
    return a, index

def vrad(z1,z2):
    '''
    Given the redshifts of the two subclusters this function calculates their
    relative line-of-sight velocity (km/s).
    Input:
    z1 = [float] redshift of subcluster 1
    z2 = [float] redshift of subcluster 2
    Output:
    v_los = [float; units:km/s] line-of-sight relative velocity of the two
        subclusters
    '''
    v1 = ((1+z1)**2-1)/((1+z1)**2+1)*c
    v2 = ((1+z2)**2-1)/((1+z2)**2+1)*c
    return numpy.abs(v1-v2)/(1-v1*v2/c**2)

def f(x,a,b):
    '''
    This is the functional form of the time since merger integral for two point
    masses.
    '''
    return 1/numpy.sqrt(a+b/x)
def TSMptpt(m_1,m_2,r_200_1,r_200_2,d_end,E):
    '''
    This function calculates the time it takes for the system to go from a
    separation of r_200_1+r_200_2 to d_end. It is based on the equations of
    motion of two point masses, which is valid in the regime where the 
    subclusters no longer overlap.
    Input:
    m_1 = [float; units: M_sun] mass of subcluster 1
    m_2 = [float; units: M_sun] mass of subcluster 2
    r_200_1 = [float; units: Mpc] NFW r_200 radius of subcluster 1
    r_200_2 = [float; units: Mpc] NFW r_200 radius of subcluster 2
    d_end = [float; units: Mpc] the final separation of interest
    E = [float; units: (km/s)^2*M_sun] the total energy (PE+KE) of the two subcluster system
    Output:
    t = [float; units: Gyr] the time it takes for the system to go from a
        separation of r_200_1+r_200_2 to d_end
    '''
    d_start = r_200_1+r_200_2
    C = G*m_1*m_2
    mu = m_1*m_2/(m_1+m_2)
    if E < 0:
        integral = scipy.integrate.quad(lambda x: f(x,E,C),d_start,d_end)[0]
        t = numpy.sqrt(mu/2)*integral/sinGyr*kminMpc        
    else:
        print 'TSMptpt: error total energy should not be > 0, exiting'
        sys.exit()
    if t < 0:
        print 'TSM < 0 encountered'    
    return t
      

def PEnfwnfw(d,m_1,rho_s_1,r_s_1,r_200_1,m_2,rho_s_2,r_s_2,r_200_2,N=100):
    '''
    This function calculates the potential energy of two truncated NFW halos.
    Input:
    d = [float; units:Mpc] center to center 3D separation of the subclusters
    m_1 = [float; units:M_sun] mass of subcluster 1 out to r_200
    rho_s_1 = [float; units:M_sun/Mpc^3] NFW scale density of subcluster 1
    r_s_1 = [float; units:Mpc] NFW scale radius of subcluster 1
    r_200_1 = [float; units:Mpc] r_200 of subcluster 1
    m_2 = [float; units:M_sun] mass of subcluster 2 out to r_200
    rho_s_2 = [float; units:M_sun/Mpc^3] NFW scale density of subcluster 2
    r_s_2 = [float; units:Mpc] NFW scale radius of subcluster 2
    r_200_2 = [float; units:Mpc] r_200 of subcluster 2
    N = [int] number of mass elements along one coordinate axis for numerical
        integration approximation
    Output:
    V_total = [float; units:(km/s)^2*M_sun] total potential energy of the two
        subcluster system
    '''
    if d >= r_200_1+r_200_2:
        V_total = -G*m_1*m_2/d
    else:
        # mass element sizes
        dr = r_200_2/N
        dt = numpy.pi/N
        i,j = numpy.meshgrid(numpy.arange(N),numpy.arange(N))
        # distance of 2nd NFW halo mass element from center of 1st NFW halo
        r = numpy.sqrt(((2*i+1)*dr/2*numpy.sin((2*j+1)*dt/2))**2+
                       (d+(2*i+1)*dr/2*numpy.cos((2*j+1)*dt/2))**2)
        #mass elements of 2nd NFW halo  
        m = 2*numpy.pi*rho_s_2*r_s_2**3*(numpy.cos(j*dt)-numpy.cos((j+1)*dt))*(1/(1+(i+1)*dr/r_s_2)-1/(1+i*dr/r_s_2)+numpy.log(((i+1)*dr+r_s_2)/(i*dr+r_s_2)))
        #determine if 2nd halo mass element is inside or outside of 1st halo
        mask_gt = r>=r_200_1
        mask_lt = r<r_200_1
        # potential energy of each mass element (km/s)^2 * M_sun
        # NFW PE
        V_nfw = numpy.sum(-4*numpy.pi*G*rho_s_1*r_s_1**3/(r[mask_lt]+dr)*(numpy.log(1+(r[mask_lt]+dr)/r_s_1)-(r[mask_lt]+dr)/(r_s_1+r_200_1))*m[mask_lt])
        # Point mass PE
        V_pt = numpy.sum(-G*m_1*m[mask_gt]/r[mask_gt])
        V_total = V_nfw+V_pt
    return V_total
    
def NFWprop(M_200,z,c):
    '''
    Determines the NFW halo related properties. Added this for the case of user
    specified concentration.
    Input:
    M_200 = [float; units:M_sun] mass of the halo. Assumes M_200 with
        respect to the critical density at the halo redshift.    
    z = [float; unitless] redshift of the halo.
    c = [float; unitless] concentration of the NFW halo.
    '''
    # CONSTANTS
    rho_cr = cosmo.rhoCrit(z)/kginMsun*minMpc**3
    #calculate the r_200 radius
    r_200 = (M_200*3/(4*numpy.pi*200*rho_cr))**(1/3.)
    del_c = 200/3.*c**3/(numpy.log(1+c)-c/(1+c))
    r_s = r_200/c
    rho_s = del_c*rho_cr   
    return del_c, r_s, r_200, c, rho_s

def MCengine(N_mc,M1,M2,Z1,Z2,D_proj,prefix,C1=None,C2=None,del_mesh=100,TSM_mesh=200):
    '''
    This is the Monte Carlo engine that draws random parameters from the
    measured distributions and then calculates the kenematic parameters of the
    merger.
    Input:
    N_mc = [int] number of successful Monte Carlo samplings to perform
    M1 = [(float, float) or (1D array of floats); units:M_sun] M_200 of 
        subcluster 1.  If (float, float) is input then it is interpreted as 
        (mu, sigma) of a Gaussian distribution. If (1D array of floats with 
        length > 2) then it is interpreted as a random sample from a parent 
        distribution.
    M2 = [(float, float) or (1D array of floats); units:M_sun] M_200 of 
        subcluster 2. Must have same format as M1.
    Z1 = [(float, float) or (1D array of floats)] redshift of subcluster 1.
        Similar input format to M1.
    Z2 = [(float, float) or (1D array of floats)] redshift of subcluster 2.
        Similar input format to M1.
    D_proj = [(float, float, float) or (1D array of floats); units:Mpc]
        projected separation of the two subclusters.  If (float, float, float)
        is input then it is interpreted as (d_proj, c1_sigma, c2_sigma), where
        d_proj is the projected distance of the two halos and c1_sigma and
        c2_sigma are the centroid errors of each halo.  If (1D array of floats
        with length > 3) then it is interpreted as a random sample from the 
        parent distribution of d_proj.
    prefix = [string] prefix to attach to all output
    C1 = [(float, float) or (1D array of floats); unitless] NFW concentration of 
        subcluster 1.  If (float, float) is input then it is interpreted as 
        (mu, sigma) of a Gaussian distribution. If (1D array of floats with 
        length > 2) then it is interpreted as a random sample from a parent 
        distribution. If None then Duffey et al. (2008) scaling relation is used
        to determine the concentration based on M1 value. Must have same format
        as M1.
    C2 = [(float, float) or (1D array of floats); unitless] NFW concentration of 
        subcluster 2. Must have same format as M2.
    del_mesh = [int] number of mass elements along one coordinate axis for
        numerical integration approximation of potential energy
    TSM_mesh = [int] number of elements along the separation axis when
        performing the numerical integration appoximation of the TSM
    Output:
    In general the 1D arrays of floats are output for each of the following
    parameters.  The Nth element of each array corresponds to the system
    properties of the Nth viable solution.
    m_1_out = [units: M_sun] M_200 of halo 1
    m_2_out = [units: M_sun] M_200 of halo 1
    z_1_out = redshift of halo 1
    z_2_out = redshift of halo 2
    alpha_out = [units: degrees] merger axis angle with respect to the plane of
        the sky. 0 corresponds to in the plane of the sky, 90 corresponds to 
        along the line-of-sight.
    d_proj_out = [units: Mpc] projected center to center separation of the
        halos in observed state
    d_3d_out = [units: Mpc] 3D center to center separation of the halos in 
        observed state
    d_max_out = [units: Mpc] 3D maximum center to center separation of the two
        halo at the apoapsis
    v_rad_obs_out = [units: km/s] relative line-of-sight velocity of the two
        halos at the observed time
    v_3d_obs_out = [units: km/s] relative 3D velocity of the two halos at the
        observed time
    v_3d_col_out = [units: km/s] relative 3D velocity of the two halos at the
        collision time
    TSM_0_out = [units: Gyr] time it took the system to reach the observed state
        from the collision state.  Assuming the system has not reached its
        apoapsis since colliding.
    TSM_1_out = [units: Gyr] time it took the system to reach the observed state
        from the collision state.  Assuming the system has passed throught the 
        apoapsis once since colliding.
    T_out = [units: Gyr] period of the system
    prob_out = This is the probability of observing the system with the randomly
       drawn paramers.  This accounts for the fact that it is more likely to
       observe the system near the apoapsis rather than with zero separation,
       due to the system spending more of its time at large separation.  Simply
       an effect of the velocity of the system.
    '''
    # Verify user input
    # Check to make sure mass inputs for each halo are the same
    N_M1 = numpy.size(M1)
    N_M2 = numpy.size(M2)
    if N_M1 != N_M2:
        print 'MCMAC.MCengine: Error, the mass inputs for the two halos must \
            be the same type and size. For example, you cannot mix input \
            format (mu, sigma) for halo1 with (1D array of floats) format \
            for halo 2. Nor can the size of the (1D array of floats) input \
            be different. This is to facilitate correct covariance handeling. \
            Exiting.' 
        sys.exit()
    
    if C1 != None or C2 != None:
        N_C1 = numpy.size(C1)
        N_C2 = numpy.size(C2)
        if N_C1 != N_M1 or N_C2 != N_M2:
            print 'MCMAC.MCengine: Error, the concentration inputs for the two \
            halos must be the same type and size as the mass inputs. For \
            example, you cannot mix input format (mu, sigma) for M1 with \
            (1D array of floats) format for C1, or vise versa. Nor can the size\
            of the (1D array of floats) input be different. This is to \
            facilitate correct covariance handeling. Exiting.'
            sys.exit()
        elif N_M1 > 2:
            print 'MCMAC.MCengine: Assuming that the order and values of the \
            C1 and C2 (1D array of floats) are correlated with the order and \
            values of the M1 and M2 (1D array of floats). This is meant to \
            maintain proper covariance.'
    
    N_D_proj = numpy.size(D_proj)
    if N_D_proj == N_M1:
        print 'MCMAC.MCengine: Assuming that the order and values of the \
        D_proj (1D array of floats) are correlated with the order and values of\
        the M1 and M2 (1D array of floats). This is meant to maintain proper\
        covariance.'
        
    N_Z1 = numpy.size(Z1)
    if N_Z1 == N_D_proj:
        print 'MCMAC.MCengine: Warning. It is currently assumed that the halo\
        redshifts and projected separation estimates are uncorrelated. \
        Covariance will not be handled correctly.'
    if N_Z1 >2 and N_Z1 == N_M1:
        print 'MCMAC,MCengine: Warning. It is currently assumed that the halo \
        redshifts and mass estimates are uncorrelated. Covariance will not be \
        handled correctly.'
        
    i = 0
    # Create the output arrays
    m_1_out = numpy.zeros(N_mc)
    m_2_out = numpy.zeros(N_mc)
    z_1_out = numpy.zeros(N_mc)
    z_2_out = numpy.zeros(N_mc)
    d_proj_out = numpy.zeros(N_mc)
    v_rad_obs_out = numpy.zeros(N_mc)
    alpha_out = numpy.zeros(N_mc)
    v_3d_obs_out = numpy.zeros(N_mc)
    d_3d_out = numpy.zeros(N_mc)
    v_3d_col_out = numpy.zeros(N_mc)
    d_max_out = numpy.zeros(N_mc)
    TSM_0_out = numpy.zeros(N_mc)
    TSM_1_out = numpy.zeros(N_mc)
    T_out = numpy.zeros(N_mc)
    prob_out = numpy.zeros(N_mc)
    
    N_d = numpy.size(D_proj)
    t_start = time.time()
    while i < N_mc:
        # Draw random mass and redshift parameters
        m_1, index = randdraw(M1)
        m_2, index = randdraw(M2,index)
        if m_1 <= 0 or m_2 <=0:
            # Discard these unphysical draws, this can happen if a Gaussian 
            # distribution is speficied since at least part of the tail will be
            # in the negative mass range.
            continue
        z_1, temp = randdraw(Z1)
        z_2, temp = randdraw(Z2)
        
        # Draw random projected separation
        if N_d == 3: #then assumed (D_proj, c1_sigma, c2_sigma)
            #draw del_r1, del_r2, phi1 and phi2, where phi is the 
            del_r_1 = D_proj[1] * numpy.random.randn()
            del_r_2 = D_proj[2] * numpy.random.randn()
            phi_1 = numpy.random.rand()*2*numpy.pi
            phi_2 = numpy.random.rand()*2*numpy.pi
            #calculate the projected separation based on drawn centers
            d_proj = numpy.sqrt((D_proj[0]-del_r_1*numpy.cos(phi_1)+del_r_2*numpy.cos(phi_2))**2+(del_r_1*numpy.sin(phi_1)+del_r_2*numpy.sin(phi_2))**2)
        elif N_d > 3: #then assume distribution array format
            if N_D_proj == N_M1:
                # assume D_proj array values and order correlated with M1 and M2
                # array values and order.
                d_proj, index = randdraw(D_proj,index)
            else:
                # assume D_proj is uncorrelated with M1 and M2
                d_proj = D_proj[numpy.random.randint(N_d)]
        elif N_d < 3:
            print 'MCMAC.MCengine: Error, D_proj parameter input array is not a valid size, exiting'
            sys.exit()
        d_proj = abs(d_proj)
        
        # Define NFW halo properties
        if C1 == None:
            # Then use scaling relation to determine concentration
            del_c_1, r_s_1, r_200_1, c_1, rho_s_1 = profiles.nfwparam_extended(m_1/1e14,z_1)
        else:
            # Then user specified concentration
            # Draw random concentration parameter
            c_1, index = randdraw(C1,index)
            del_c_1, r_s_1, r_200_1, c_1, rho_s_1 = NFWprop(m_1,z_1,c_1)
        if C2 == None:
            del_c_2, r_s_2, r_200_2, c_2, rho_s_2 = profiles.nfwparam_extended(m_2/1e14,z_2)
        else:
            c_2, index = randdraw(C2,index)
            del_c_2, r_s_2, r_200_2, c_2, rho_s_2 = NFWprop(m_2,z_2,c_2)
        
        # Reduced mass
        mu = m_1*m_2/(m_1+m_2)        
        
        # Calculate potential energy at time of collision
        PE_col = PEnfwnfw(0,m_1,rho_s_1,r_s_1,r_200_1,m_2,rho_s_2,r_s_2,r_200_2,N=del_mesh)
        
        # Calculate maximum 3D free-fall velocity
        # Note that this assumes that the clusters begin with infinite
        # separation and their observed mass. It also guarntees that all
        # solutions will be for a bound system.
        v_3d_max = numpy.sqrt(-2/mu*PE_col)
        
        # Calculate observed radial velocity 
        v_rad_obs = vrad(z_1,z_2)
        
        # Draw a random alpha, merger axis angles with respect to the sky
        ALPHA = numpy.random.rand()*numpy.pi/2
        #Since not all alpha are equally likely.
        alpha = (1-numpy.cos(ALPHA))*numpy.pi/2
        
        # Calculate the 3D velocity at observed time
        v_3d_obs = v_rad_obs/numpy.sin(alpha)
        
        if v_3d_obs > v_3d_max:
            # then the randomly chosen alpha is not valid
            continue
        
        # Calculate the 3D separation
        d_3d = d_proj/numpy.cos(alpha)
        
        # Calculate the potential energy at observed time
        PE_obs = PEnfwnfw(d_3d,m_1,rho_s_1,r_s_1,r_200_1,m_2,rho_s_2,r_s_2,r_200_2,N=del_mesh)
        
        # Calculate the 3D velocity at collision time
        v_3d_col = numpy.sqrt(v_3d_obs**2+2/mu*(PE_obs-PE_col))
        
        if v_3d_col > v_3d_max:
            # then the randomly chosen alpha is not valid
            continue
        
        # Total Energy
        E = PE_obs+mu/2*v_3d_obs**2
        
        # Calculate PE from d = 0 to r_200_1+r_200_2
        del_TSM_mesh = (r_200_1+r_200_2)/(TSM_mesh-1)
        d = numpy.arange(0.00001,(r_200_1+r_200_2)+del_TSM_mesh,del_TSM_mesh)
        PE_array = numpy.zeros(TSM_mesh)
        for j in numpy.arange(TSM_mesh):
            PE_array[j] = PEnfwnfw(d[j],m_1,rho_s_1,r_s_1,r_200_1,m_2,rho_s_2,r_s_2,r_200_2,N=del_mesh)        

        # Calculate d_max        
        if E >= -G*m_1*m_2/(r_200_1+r_200_2):
            # then d_max > r_200_1+r_200_2
            d_max = -G*m_1*m_2/E
        else:
            # d_max <= r_200_1+r_200_2
            ## find closet d value less than d_max
            #idx=(numpy.abs(PE_array-E)).argmin()
            #d_max = d[idx]
            # determine d_max, ensuring that E is always > PE
            mask_tmp = E >= PE_array
            d_max = d[mask_tmp][-1]
            
        # Calculate TSM_0
        if d_3d >= r_200_1+r_200_2:
            # then halos no longer overlap
            # calculate the time it takes to go from d=0 to r_200_1+r_200_2
            TSM_0a = numpy.sum(del_TSM_mesh/numpy.sqrt(2/mu*(E-PE_array))*kminMpc/sinGyr)
            # calculate th time it takes to go from d = r_200_1+r_200_2 tp d_3d
            TSM_0b = TSMptpt(m_1,m_2,r_200_1,r_200_2,d_3d,E)
            TSM_0 = TSM_0a+TSM_0b
        else:
            # then d_3d < r_200_1+r_200_2, halos always overlap
            # calculate the time it takes to go from d=0 to d_3d            
            mask = d <= d_3d
            TSM_0 = numpy.sum(del_TSM_mesh/numpy.sqrt(2/mu*(E-PE_array[mask]))*kminMpc/sinGyr)
        
        # Check that TSM_0 < Age of Universe at (z_1+z_2)/2
        age = cosmo.age((z_1+z_2)/2)
        if TSM_0 > age:
            # unlikely that this system could occur
            continue
        
        # Calculate period
        if E >= -G*m_1*m_2/(r_200_1+r_200_2):
            # then d_max > r_200_1+r_200_2
            # calculate th time it takes to go from d = r_200_1+r_200_2 to d_max
            if d_3d < r_200_1+r_200_2:
                #then TSM_0a has not been previously defined
                TSM_0a = numpy.sum(del_TSM_mesh/numpy.sqrt(2/mu*(E-PE_array))*kminMpc/sinGyr)
            TSM_0b = TSMptpt(m_1,m_2,r_200_1,r_200_2,d_max,E)
            T = 2*(TSM_0a+TSM_0b)
        else:
            # then d_max < r_200_1+r_200_2
            # calculate the time it takes to go from d=0 to d_max           
            mask = d <= d_max
            T = 2*numpy.sum(del_TSM_mesh/numpy.sqrt(2/mu*(E-PE_array[mask]))*kminMpc/sinGyr)        
            
        # Calculate probability of d_3d
        prob = TSM_0/(T/2)
        
        # Calculate TSM_1
        TSM_1 = T-TSM_0        

        if TSM_0 < 0:
            print 'TSM < 0 encountered'
        
        # Write calculated merger parameters
        m_1_out[i] = m_1
        m_2_out[i] = m_2
        z_1_out[i] = z_1
        z_2_out[i] = z_2
        d_proj_out[i] = d_proj
        v_rad_obs_out[i] = v_rad_obs
        alpha_out[i] = alpha*180/numpy.pi
        v_3d_obs_out[i] = v_3d_obs
        d_3d_out[i] = d_3d
        v_3d_col_out[i] = v_3d_col
        d_max_out[i] = d_max
        TSM_0_out[i] = TSM_0
        TSM_1_out[i] = TSM_1
        T_out[i] = T
        prob_out[i] = prob        
        
        i+=1
        
        # Estimate calculation time
        if i== 10 or i == 100 or i%1000 == 0:
            del_t = time.time()-t_start
            t_total = N_mc*del_t/i
            eta = (t_total-del_t)/60
            print 'Completed Monte Carlo iteration {0} of {1}.'.format(i,N_mc)
            print '~{0:0.0f} minutes remaining'.format(eta)

    # Pickle the results of the MC analysis
    filename = prefix+'_m_1.pickle'
    F = open(filename,'w')
    pickle.dump(m_1_out,F)
    F.close()
    filename = prefix+'_m_2.pickle'
    F = open(filename,'w')
    pickle.dump(m_2_out,F)
    F.close()
    filename = prefix+'_z_1.pickle'
    F = open(filename,'w')
    pickle.dump(z_1_out,F)
    F.close()    
    filename = prefix+'_z_2.pickle'
    F = open(filename,'w')
    pickle.dump(z_2_out,F)
    F.close()
    filename = prefix+'_d_proj.pickle'
    F = open(filename,'w')
    pickle.dump(d_proj_out,F)
    F.close()
    filename = prefix+'_v_rad_obs.pickle'
    F = open(filename,'w')
    pickle.dump(v_rad_obs_out,F)
    F.close()
    filename = prefix+'_alpha.pickle'
    F = open(filename,'w')
    pickle.dump(alpha_out,F)
    F.close()
    filename = prefix+'_v_3d_obs.pickle'
    F = open(filename,'w')
    pickle.dump(v_3d_obs_out,F)
    F.close()  
    filename = prefix+'_d_3d.pickle'
    F = open(filename,'w')
    pickle.dump(d_3d_out,F)
    F.close()    
    filename = prefix+'_v_3d_col.pickle'
    F = open(filename,'w')
    pickle.dump(v_3d_col_out,F)
    F.close()    
    filename = prefix+'_d_max.pickle'
    F = open(filename,'w')
    pickle.dump(d_max_out,F)
    F.close() 
    filename = prefix+'_TSM_0.pickle'
    F = open(filename,'w')
    pickle.dump(TSM_0_out,F)
    F.close()  
    filename = prefix+'_TSM_1.pickle'
    F = open(filename,'w')
    pickle.dump(TSM_1_out,F)
    F.close()    
    filename = prefix+'_T.pickle'
    F = open(filename,'w')
    pickle.dump(T_out,F)
    F.close()  
    filename = prefix+'_prob.pickle'
    F = open(filename,'w')
    pickle.dump(prob_out,F)
    F.close()    
    
    
    return m_1_out,m_2_out,z_1_out,z_2_out,d_proj_out,v_rad_obs_out,alpha_out,v_3d_obs_out,d_3d_out,v_3d_col_out,d_max_out,TSM_0_out,TSM_1_out,T_out,prob_out