plotCorrelationMap11.py

#!/usr/bin/env python
"""
By David French (frenchd@astro.wisc.edu)

$Id: plotCorrelationMap11.py, v 11.3 10/03/16

This program takes in a list of AGN targets and generates an environment map (i.e. nearby
galaxies) for each. 


Adapted from:
'
    plotAGNCorrelationSquare9.py, v 9.0 02/12/2014

    Make a square plot showing the positions of correlated galaxies around a target AGN

    Starts with the older 'plotAGNCorrelationSquare6.py' and makes adjustments for the current
    data.

    v8: include inclinations in galaxy annotation tags

    v9: plots all targets in the input file, and does everything better faster
'

v1: Removes group vs field galaxy bits

v2: Draw ellipses to denote inclination and position angle at the same time

v3: Make nicer plots for poster (09/10/14)

v4: change the labeling for individual plots (09/11/14)

v5: adopt to make plots for general correlation tables

v6: include option to make skyplot as well (plot just the ra and dec of all objects)
    This was final AAS poster version, made targetmaps3
    
v7: More general version - 06/19/15

v7.1: Now saves into 'ambiguous' or 'associated' folders depending on the 'include' answer
    - 06/25/15

7.2: Reverse PA direction. Now PA is measured clockwise from N (up) to E (on right). This
    changes angle = float(PA[i]) to angle = float(-PA[i]) in plotting
    
    Also changed the name annotation slightly, to xy=(-60, size[i]/6) from 
    xy=(-65, size[i]/4), and unknown type marker size increased to size[i]*4
    
    E type galaxies now have no edge, spiral have black edges, unknown are diamonds
    
    - 07/07/15
    
    - unfortunately, this version was still named "plotCorrelationMap7_1.py".
    

7.3: Sort the map table by the quantity: abs(dif_velocity)*impact parameter/virial radius
    - low means the galaxy in question close in both velocity and distance. These get
     listed first.

7.4: Same as above, but tweaking the "likelihood" parameter to the following:
    abs(dif_velocity) * (impact parameter - virial radius)

7.5: Now the likelihood parameter is:
      abs(dif_velocity) * (impact parameter - virial radius) / (virial radius)**2
      
      if this value is >2 times the next largest, call it associated and put it in that folder
      
      11/6/15: (approximate date) Now use a gaussian form likelihood parameter, and take
      the velocity difference into account as well (normalized to 200 km/s)
        -looks like: exp((impact/virial)^2) * exp((vel_dif/200)^2) - Targetmaps12
      
      11/11/15: continued with more tweaks to try to get satellite galaxies to be rated
            lower. Targetmaps13 uses diam^1.4 instead of a virial radius
            
            Targetmaps14 uses diam^1.5
            
    11/15/15: Targetmaps15 uses diam^1.5, and multiplies likelihood by 2 if the virial
            radius is larger than the impact parameter
            
            Targetmaps16 requires 3* greater likelihood for associated tag
            
            Targetmaps17 back to standard virial radius, but with the above tweaks for 
            multiplying likelihood *2 and requiring 3* difference
            
    11/17/15: Targetmaps18 back to diam^1.5, and now with 5* difference required
    
    11/19/15: Targetmaps19 uses standard r_vir, but requires a 5* difference
    
    
v8.0: Computes both standard R_vir and diam^1.5 results. Adds the results automatically to
LG_correlation_combined5.csv, 

v9.0: velocity colormap in units of delta-v instead of absolute velocity units

v9.1: 'include' now requires a hard L limit (L>=0.001 etc)

v10: correlate with galaxy v_helio, not vcorr - using correlateSingle7 now
    - (03/24/16 - 4/13/16)
    
v11: moved to .../inclination/git_inclination/ so it can be updated with git, and 
    now including rc tex formating for the plots (05/12/2016)
    
v11.1: updates for the newest round of sightlines in LG_correlation_combined5_10.csv
       - Made LG_correlation_combined5_11.csv and targetmaps34 (07/06/16)
       
v11.2: minor formatting updates. (8/08/16)

v11.3: more minor updates - make the tick labels NOT bold (10/03/16)
    - make LG_correlation_combined5_11_25cut_edit5.csv and targetmaps37/

v11.3: more minor updates - fix the colorbar bolded ticks problem (10/10/16)
    - remake LG_correlation_combined5_11_25cut_edit5.csv and targetmaps37/
"""

import sys
import os
import csv
import string
import math
import ast
from pylab import *
import correlateSingle8 as correlateSingle
from matplotlib.patches import Ellipse
# from astropy.coordinates import SkyCoord
# from astropy import units as u
from utilities import *
import getpass

from matplotlib import rc
fontScale = 18
rc('text', usetex=True)
rc('font', size=18,family='serif',weight='medium')
rc('xtick.major',size=8,width=0.6)
rc('xtick.minor',size=5,width=0.6)
rc('ytick.major',size=8,width=0.6)
rc('ytick.minor',size=5,width=0.6)
rc('xtick',labelsize = fontScale)
rc('ytick',labelsize = fontScale)
rc('axes',labelsize = fontScale)
rc('xtick', labelsize = fontScale)
rc('ytick',labelsize = fontScale)
# rc('font', weight = 450)
# rc('axes',labelweight = 'bold')
rc('axes',linewidth = 1,labelweight='normal')

# from matplotlib.patches import Ellipse

    
################################################################


    
def buildTargetList(file,AGNheader,velocityHeader):
    # builds a list consisting of tuples: (AGNname, velocity, include)
    #
    # requires a column called 'galaxyName' in the file. Only includes targets for which
    # 'galaxyName' != 'x'
    #
    # - AGN name is the name of the target
    # - velocity is the velocity of the absorption feature centroid in km/s
    # - include is a boolean indicating whether the results should be saved into the
    #   'ambiguous' folder (False), or the 'associated' folder (True)
    
    f = open(file,'rU')
    reader = csv.DictReader(f)
    targetList = []
    for l in reader:
        AGNname = l[AGNheader]
        include = l['include']
        if include == 'yes':
            include = True
        else:
            include = False
        
        # velocityHeader is the name of the columns containing the center velocity, or 
        # it can also be just a number designating the center velocity
        if isNumber(velocityHeader):
            vel = velocityHeader
        else:
            vel = l[velocityHeader]
            
        galaxyName = l['galaxyName']
        if isNumber(vel) and galaxyName !='x':
            pair = (AGNname,int(vel), include)
            targetList.append(pair)
    
    f.close()
    return targetList


def buildFullTargetList(file,AGNheader,velocityHeader):
    # makes a targetmap for every absorption line, regardless of whether there is an
    # associated galaxy name or not, so long as the line velocity entry is a number
    #
    # builds a list consisting of tuples: (AGNname, velocity, include)
    #
    # - AGN name is the name of the target
    # - velocity is the velocity of the absorption feature centroid in km/s
    # - include is a boolean indicating whether the results should be saved into the
    #   'ambiguous' folder (False), or the 'associated' folder (True)
    
    f = open(file,'rU')
    reader = csv.DictReader(f)
    targetList = []
    for l in reader:
        AGNname = l[AGNheader]
        include = l['include']
        if include == 'yes':
            include = True
        else:
            include = False
        
        # velocityHeader is the name of the columns containing the center velocity, or 
        # it can also be just a number designating the center velocity
        if isNumber(velocityHeader):
            vel = velocityHeader
        else:
            vel = l[velocityHeader]
            
        if isNumber(vel):
            pair = (AGNname,int(vel),include)
            targetList.append(pair)
    
    f.close()
    return targetList
        
        
def main():
    # main function to create targetmaps around selected sightlines
        
    counter = 0
    AGNList = []
    rowList = []
    masterVirList = []
    masterCustomList = []
    
    # max impact parameter to use
    maxSep = 500
    
    # +/- galaxy velocity to search within around absorption velocity
    velocityWindow = 400
    
    # minimum galaxy velocity to include (False to ignore)
    minVcorr = 500
    
    # minimum galaxy size to include (False to ignore)
    minSize = False
    
    # minimum separation in km/s between the redshift of the AGN and the galaxy (False to ignore) 
    agnSeparation = False
    
    # include name tags on galaxies? They don't scale very well...
    includeNameTags = True
    
    # include a title on the plots?
    includeTitle = False
    
    # also make a plot of just real positions of galaxies on the sky in RA and Dec coords?
    includeSkyPlot = False
    
    # Save the map plots?
    saveMaps = True
    
    # Save the individual map plot tables?
    saveMapTables = True
    
    # Save the full results with "include" tags?
    saveResults = False
    
    # 2nd place galaxy likelihood * rigor <= 1st place galaxy for 'include'
    rigor = 5
    
    # hard limit for likelihood
    l_min = 0.01
    
    # bypass l_min for lone galaxies? (i.e. include lone galaxies no matter what likelihood is)
    loner = False
    
    # sort results into /associated/, ~/ambiguous/, and ~/nonassociated/ folders?
    # if True, these folders must already exist
    # if False, puts all the files into saveDirectory as set below
    sortIntoFolders = True
    
    
    # where to save figures and tables
    user = getpass.getuser()
    if user == "David":
        targetFile = '/Users/David/Research_Documents/inclination/git_inclination/LG_correlation_combined5_11_25cut_edit4.csv'
        saveDirectory = '/Users/David/Research_Documents/inclination/git_inclination/targetmaps38/'
        outputFile = '/Users/David/Research_Documents/inclination/git_inclination/LG_correlation_combined5_11_25cut_edit5.csv'

    elif user == "frenchd":
        targetFile = '/usr/users/frenchd/inclination/git_inclination/LG_correlation_combined5_11_25cut_edit4.csv'
#         saveDirectory = '/usr/users/frenchd/inclination/git_inclination/targetmaps38/'
        saveDirectory = '/usr/users/frenchd/iraf/NGC5364/'
#         outputFile = '/usr/users/frenchd/inclination/git_inclination/LG_correlation_combined5_11_25cut_edit5.csv'
        outputFile = '/usr/users/frenchd/iraf/NGC5364_correlationMap.csv'

    else:
        print "Unknown user: ",user
        sys.exit()
    
    # what are the column names in this file for the AGN name and absorption velocity?
    AGNheader = 'AGNname'
    velocityHeader = 'Lya_v'

    # targets from a file, use this:
#     targets = buildFullTargetList(targetFile,AGNheader,velocityHeader)

    
    # or build up a custom list of AGN names and absorption velocities here:
    targets = [('SDSSJ135726.27+043541.4',1094.0,True)]

    
    c = 0
    for i in targets:
        # find AGN environment using the imported version of correlateSingle
        AGNname,center,include = i
        correlation = correlateSingle.correlateTarget(AGNname, maxSep, agnSeparation, minVcorr, minSize, slow=False)
        galaxyInfo = correlation[AGNname]
        
        print '{0} = {1}'.format(AGNname,len(galaxyInfo))
        print
        
#         galaxyInfo.sort()

        # instantiate some lists for later
        galaxyNames = []
        separations = []
        positions = []
        plotPositionsRA = []
        plotPositionsDec = []
        plotRA = []
        plotDec = []
        plotAGNposition = []
        plotSizes = []
        plotVelocity = []
        PA = []
        inc = []
        typeList = []
        infoDict = {}
        
        # loop through the returned galaxy environment data, making calculations and
        # populating lists as we go
        counter = 0
        for row in galaxyInfo:
            counter+=1
            vhel, galaxyRow = row
            AGNposition = eval(str(galaxyRow['degreesJ2000RA_DecAGN']))

            
            # crop off results that fall out of the 'velocityWindow' parameter
            # 'velocityWindow' = a cut in velocity space 
            if counter <=5000 and float(vhel)-velocityWindow <= center and float(vhel)+velocityWindow >= center:
                separation = galaxyRow['impactParameter (kpc)']
                galaxyName = galaxyRow['galaxyName']
                galaxyPosition = eval(str(galaxyRow['degreesJ2000RA_DecGalaxy']))
                AGNposition = eval(str(galaxyRow['degreesJ2000RA_DecAGN']))
                galaxyDist = galaxyRow['distGalaxy (Mpc)']
                group = eval(str(galaxyRow['groups_dist_std (Mpc)']))
                major,minor = eval(str(galaxyRow['linDiameters (kpc)']))
                galaxyVcorr = galaxyRow['vcorrGalaxy (km/s)']
                galaxyVel = galaxyRow['radialVelocity (km/s)']
                morphology = galaxyRow['morphology']
                RC3Type = galaxyRow['RC3type']
                RC3PA = galaxyRow['RC3pa (deg)']
                RC3inc = galaxyRow['RC3inc (deg)']
                positionAngle = galaxyRow['positionAngle (deg)']
                inclination = galaxyRow['inclination (deg)']
                azimuth = galaxyRow['azimuth (deg)']
#                 include = galaxyRow['include']
                
                if isNumber(inclination):
                    inclination = round(eval(inclination),0)
                elif isNumber(RC3inc):
                    inclination = round(eval(RC3inc),0)

                if isNumber(galaxyDist):
                    positions.append(galaxyPosition)
                    separations.append(float(separation))
                    
                    if isNumber(positionAngle):
                        pa = float(positionAngle)
                    elif isNumber(RC3PA):
                        pa = float(RC3PA)
                    else:
                        pa = 0
                        
                    #find plot placement
                    # calculate angular separations in ra and dec to determine positions on chart w.r.t. target AGN
                    gRA,gDec = float(galaxyPosition[0]),float(galaxyPosition[1])
                    agnRA,agnDec = float(AGNposition[0]),float(AGNposition[1])
                    
                    # calculate separation in RA only
#                     dRA = correlateSingle.calculateImpactParameter_slow(gRA,agnDec,agnRA,agnDec,galaxyDist)
                    dRA = correlateSingle.calculateImpactParameter(gRA,agnDec,agnRA,agnDec,galaxyDist)

                    
                    # calculate separation in Dec only
#                     dDec = correlateSingle.calculateImpactParameter_slow(agnRA,gDec,agnRA,agnDec,galaxyDist)
                    dDec = correlateSingle.calculateImpactParameter(agnRA,gDec,agnRA,agnDec,galaxyDist)

                    
                    # add signs back into physical impact parameters
                    if gRA < agnRA:
                        dRA = -dRA
                        
                    # 'edge' effects
                    if gRA >= 359.0 and agnRA <= 1.0:
                        dRA = -dRA
                        
                    if gDec < agnDec:
                        dDec = -dDec
                        
                        
                    # calculate size by finding radius
                    noSize = False
                    if isNumber(major):
                        averageSize = float(major)
                    elif isNumber(minor) and isNumber(inclination):
                        averageSize = float(minor) / math.cos(float(inclination) * math.pi/180)
                    elif isNumber(minor) and not isNumber(inclination) and not isNumber(major):
                        averageSize = float(minor)
                        inclination = 0
                    else:
                        averageSize = 2
                        inclination = 0
                        noSize = True
                                            
                    if not isNumber(galaxyVcorr):
                        galaxyVcorr = 0
                        galaxyVel = 0

                    localType = 'x'
                    rc3Type = RC3Type.lower()
                    r = rc3Type[:3]
                    morphology = morphology.lower()
                    m = morphology[:3]
                    if morphology != 'x':
                        if bfind(m,'s'):
                            if not bfind(m,'s0'):
                                # straight spiral type
                                localType = 's'
                            else:
                                # lenticular or S0 type
                                localType = 'e'
                        elif bfind(m,'e') or bfind(m,'dwarf') or bfind(m,'pec'):
                            localType = 'e'
                        elif bfind(m,'len'):
                            localType = 'e'
                        
                        elif bfind(m,'Ir') or bfind(m,'Im') or bfind(m,'I ') or bfind(m,'IA'):
                            localType = 's'
                        
                        else:
                            localType = 'x'
                            
                    # try RC3 types
                    elif rc3Type !='x':
                        if bfind(r,'s'):
                            if not bfind(r,'s0'):
                                # straight spiral type
                                localType = 's'
                            else:
                                # lenticular or S0 type
                                localType = 'e'
                        elif bfind(r,'e') or bfind(r,'dwarf') or bfind(r,'pec'):
                            localType = 'e'
                        elif bfind(r,'len'):
                            localType = 'e'
                            
                        elif bfind(r,'Ir') or bfind(r,'Im') or bfind(r,'I ') or bfind(r,'IA'):
                            # irregular type
                            localType = 's'
                        else:
                            localType = 'x'
                            
                    # Latex format galaxyName
                    galaxyName = r'{0}'.format(galaxyName)
                    galaxyName = galaxyName.replace('_','\_')
                            
                    infoDict[galaxyName]=galaxyRow
                    
                    galaxyType = {'s':2,'e':1,'x':3}

                    
                    if noSize:
                        galaxyNames.append('*'+galaxyName)
                    else:
                        galaxyNames.append(galaxyName)
                    
                        
                    plotPositionsRA.append(dRA)
                    plotPositionsDec.append(dDec)
                    plotSizes.append(averageSize)
                    plotVelocity.append(float(galaxyVel))
                    typeList.append(galaxyType[localType])
                    PA.append(pa)
                    inc.append(inclination)
                    plotRA.append(gRA)
                    plotDec.append(gDec)
                    plotAGNposition.append(AGNposition)
                    

        if len(galaxyNames) >=1:                        
            ##################################
            print 'starting first plot...'
            # color map: 
            colmap = cm.RdBu
            
            x = arange(len(galaxyNames))+1
            fig = figure(figsize=(9,7))
            ax = fig.add_subplot(111)
            width = 0.30
            
            # format the axes:
            #
            # x-axis
            majorLocator   = MultipleLocator(100)
            majorFormatter = FormatStrFormatter(r'$\rm %d$')
            minorLocator   = MultipleLocator(50)
            ax.xaxis.set_major_locator(majorLocator)
            ax.xaxis.set_major_formatter(majorFormatter)
            ax.xaxis.set_minor_locator(minorLocator)
        
            # y axis
            majorLocator   = MultipleLocator(100)
            majorFormatter = FormatStrFormatter(r'$\rm %d$')
            minorLocator   = MultipleLocator(50)
            ax.yaxis.set_major_locator(majorLocator)
            ax.yaxis.set_major_formatter(majorFormatter)
            ax.yaxis.set_minor_locator(minorLocator)            
            
#             ax.xaxis.set_tick_params(labelweight='normal')
#             ax.yaxis.set_tick_params(labelweight='normal')
            
            # scale sizes and velocities
            maxSize = 300
            minSize = 80
            largest = float(max(plotSizes))
            smallest = float(min(plotSizes))
            
            newSizes = []
        
            # multiply sizes of galaxies by 10
            for s in plotSizes:
                new = s*10
                newSizes.append(new)
            
            vmaxVal = 400
            vminVal = -400

            # +/- 400 km/s around the center
            largestVelocity = velocityWindow
            smallestVelocity = -velocityWindow

            newVelocities = []
            # check if there's more than one
            for v in plotVelocity:
                # convert to delta-v = v_absorber - v_galaxy (neg = absorber is blueward of galaxy)
                velocity = center - v
#                 newVelocity = ((float(velocity) - smallestVelocity)/(largestVelocity-smallestVelocity)) * (vmaxVal-0)+0
#                 newVelocities.append(newVelocity)

                newVelocities.append(velocity)

            
#             rounding = -1
#             step = round(((largestVelocity-smallestVelocity)/vmaxVal),rounding)
#             ticks = arange(int(round(smallestVelocity,-2)),int(round(largestVelocity,-2))+int(step),int(step))

            rounding = -1
            step = 50
            ticks = arange(-velocityWindow,velocityWindow+step,int(step))
            print 'ticks: ',ticks
            
            norm = matplotlib.colors.Normalize(vmin = vminVal, vmax = vmaxVal)
            m = matplotlib.cm.ScalarMappable(norm=norm, cmap=colmap)

            numberGalaxies = len(newVelocities)
            galaxyIndices = range(numberGalaxies)
            
            # make ellipses to indicate position angle and inclination at once
            if len(plotVelocity) !=0:
                mapping = dict((i,n) for i,n in zip(galaxyIndices,galaxyNames))
                
                for i in range(numberGalaxies):
                    if typeList[i] == 1 or typeList[i] == 2:
                        # typeList[i] == 2 is a spiral, 1 is E type
                        print 'i: ',i,galaxyIndices[i]
                        e = Ellipse(xy=(plotPositionsRA[i],plotPositionsDec[i]),\
                        width=newSizes[i] * math.cos(float(inc[i]) * math.pi/180)/2, height=float(newSizes[i])/2, angle=float(-PA[i]))
                        
                        # no transparency
                        e.set_alpha(0.9)
                        
                        if typeList[i] == 2:
                            # spiral - edge color is black
                            
                            ax.add_artist(e)
                            e.set_facecolor(m.to_rgba(newVelocities[i]))
                            e.set_edgecolor('black')
                        
                        if typeList[i] == 1:
                            # elliptical - edge color is same as galaxy
                            
                            ax.add_artist(e)
                            e.set_facecolor(m.to_rgba(newVelocities[i]))
                            e.set_edgecolor(m.to_rgba(newVelocities[i]))
                        
                    else:
                        ax.scatter(plotPositionsRA[i],plotPositionsDec[i],s=newSizes[i]*4,c=newVelocities[i],vmin=vminVal,vmax=vmaxVal,marker=(2,1,float(-PA[i])),lw=1,cmap=colmap)
                        
                    # annotate with galaxy names if includeNameTags == True
                    if includeNameTags:
                        if galaxyNames[i].find('*')!=-1:
                            # this indicates no size data is available
                            plt.annotate('*'+str(galaxyNames[i]),xy=(plotPositionsRA[i],plotPositionsDec[i]),xytext=(-60,newSizes[i]/6),textcoords='offset points',size=8)

                        else:
                            plt.annotate(galaxyNames[i],xy=(plotPositionsRA[i],plotPositionsDec[i]),xytext=(-60,newSizes[i]/6),textcoords='offset points',size=8)
            
            plot1 = ax.scatter(plotPositionsRA,plotPositionsDec,s=0,c=newVelocities,vmin=vminVal,vmax=vmaxVal,marker='.',cmap=colmap)

#########################################################################################
#########################################################################################
            # prepare all the stuff for writing. Deciding to write happens later

            fieldnames = ('AGNname',\
            'center',\
            'galaxyName',\
            'environment',\
            'degreesJ2000RA_DecAGN',\
            'degreesJ2000RA_DecGalaxy',\
            'likelihood',\
            'likelihood_1.5',\
            'virialRadius',\
            'd^1.5',\
            'impactParameter (kpc)',\
            'redshiftDistances',\
            'vcorrGalaxy (km/s)',\
            'radialVelocity (km/s)',\
            'vel_diff',\
            'distGalaxy (Mpc)',\
            'AGN S/N',\
            'majorAxis (kpc)',\
            'minorAxis (kpc)',\
            'inclination (deg)',\
            'positionAngle (deg)',\
            'azimuth (deg)',\
            'RC3flag',\
            'RC3type',\
            'RC3inc (deg)',\
            'RC3pa (deg)',\
            'morphology',\
            'final_morphology',\
            'galaxyRedshift')
            
            virList = []
            customList = []
            environment = numberGalaxies
            
            for number in range(numberGalaxies):
                if mapping[number].find("*")!=-1:
                    infoRow = infoDict[mapping[number].strip('*')]
                else:
                    infoRow = infoDict[mapping[number]]
                    
                inc = infoRow['inclination (deg)']
                RC3inc = infoRow['RC3inc (deg)']
                if isNumber(inc):
                    inc = round(float(inc),1)
                if isNumber(RC3inc) and str(RC3inc) != '-99':
                    RC3inc = round(float(RC3inc),1)
                    
                pa = infoRow['positionAngle (deg)']
                RC3pa = infoRow['RC3pa (deg)']
                if isNumber(pa):
                    pa = round(float(pa),1)
                if isNumber(RC3pa) and str(RC3pa) != '-99':
                    RC3pa = round(float(RC3pa),1)
                    
                az = infoRow['azimuth (deg)']
                if isNumber(az):
                    az = round(float(az),1)
                else:
                    az = 'x'
                    
                vcorr = infoRow['vcorrGalaxy (km/s)']
                vhel = infoRow['radialVelocity (km/s)']
                if isNumber(vcorr):
                    vcorr = round(float(vcorr),1)
                    vhel = round(float(vhel),1)
                    vel_dif = vhel - float(center)
                else:
                    vcorr = 'x'
                    vhel = 'x'
                    vel_dif = 'x'
                
                major,minor = eval(infoRow['linDiameters (kpc)'])
                if isNumber(major):
                    major = round(float(major),1)
                if isNumber(minor):
                    minor = round(float(minor),1)
    
                agn_sn = 'x'
                corrected_az = az
                redshiftDistance = 'x'
                impact = float(infoRow['impactParameter (kpc)'])
                
                if isNumber(major):
                    rVir = calculateVirialRadius(major)

                    # try this "sphere of influence" value instead
                    m15 = major**1.5

                    # first for the virial radius
                    likelihood = math.exp(-(impact/rVir)**2) * math.exp(-(vel_dif/200.)**2)
                    
                    if rVir>= impact:
                        likelihood = likelihood*2
                        
                    # then for the second 'virial like' m15 radius
                    likelihoodm15 = math.exp(-(impact/m15)**2) * math.exp(-(vel_dif/200.)**2)
                    
                    if m15>= impact:
                        likelihoodm15 = likelihoodm15*2
                        
                    # should be like 33% at 1R_v, and linear down after that.
                    # 66% at 0.5R_v, something quadratic, or logarithmic
                    # look up the "sphere of influence" for where the probability drops
                    # down to 10%
                    
                    # use M/L, surpress environment by M/L ratios
                    # also include velocity difference to make the virial radius 3D - 
                    # really a 3D impact parameter
                    
                else:
                    likelihood = -99
                    likelihoodm15 = -99
                    rVir = -99
                    m15 = -99
                
                objectInfoList = [AGNname,\
                center,\
                infoRow['galaxyName'],\
                environment,\
                infoRow['degreesJ2000RA_DecAGN'],\
                infoRow['degreesJ2000RA_DecGalaxy'],\
                likelihood,\
                likelihoodm15,\
                rVir,\
                m15,\
                impact,\
                redshiftDistance,\
                vcorr,\
                vhel,\
                vel_dif,\
                infoRow['distGalaxy (Mpc)'],\
                agn_sn,\
                major,\
                minor,\
                inc,\
                pa,\
                az,\
                infoRow['RC3flag'],\
                infoRow['RC3type'],\
                RC3inc,\
                RC3pa,\
                infoRow['morphology'],\
                infoRow['morphology'],\
                infoRow['galaxyRedshift']]
                
                virList.append([likelihood,objectInfoList])
                customList.append([likelihoodm15,objectInfoList])
                
            # now sort the list of galaxies by likelihood:
            sorted_virList = sorted(virList,reverse=True)
            sorted_cusList = sorted(customList,reverse=True)
            
            print AGNname,' - ',center,': ', sorted_virList
            print 'and : ', sorted_cusList
            
            # if there are enough galaxies, compare the best two : VIRIAL
            if len(sorted_virList) >=2:
                first_vir = sorted_virList[0][0]
                second_vir = sorted_virList[1][0]
            
                # the most likely galaxy must have rigor * the likelihood of the #2 galaxy
                if second_vir*rigor <= first_vir:
                    # larger than the min?
                    if first_vir >= l_min:
                        include_vir = True
                    else:
                        include_vir = False
                else:
                    include_vir = False
            
            
            # include lone galaxies if loner = TRUE, otherwise same l_min requirement
            elif len(sorted_virList) == 1:
                first_vir = sorted_virList[0][0]
                
                if loner:
                    # include by default
                    include_vir = True
                    
                elif first_vir >= l_min:
                    # include only if above l_min threshold
                    print 'first_vir >= l_min: ',first_vir
                    print 'AGNname, center = ',AGNname, center
                    include_vir = True
                else:
                    include_vir = False
                    print 'first_vir < l_min: ',first_vir
                    print 'AGNname, center = ',AGNname, center
            
            else:
                # if there are no galaxies, don't include
                include_vir = False
                

            # if there are enough galaxies, compare the best two : CUSTOM
            if len(sorted_cusList) >=2:
                first_cus = sorted_cusList[0][0]
                second_cus = sorted_cusList[1][0]
            
                # the most likely galaxy must have rigor * the likelihood of the #2 galaxy
                if second_cus*rigor <= first_cus:
                    # above the hard lower limit?
                    if first_cus >= l_min:
                        include_cus = True
                    else:
                        include_cus = False
                else:
                    include_cus = False                
            
            # include lone galaxies if loner = TRUE
            elif len(sorted_cusList) == 1:
                first_cus = sorted_cusList[0][0]
                
                if loner:
                    # include by default
                    include_cus = True
                    
                elif first_cus >= l_min:
                    # include only if above l_min threshold
                    print 'first_cus >= l_min: ',first_cus
                    print 'AGNname, center = ',AGNname, center
                    include_cus = True
                else:
                    # don't include
                    include_cus = False
                    print 'first_cus < l_min: ',first_cus
                    print 'AGNname, center = ',AGNname, center
            
            else:
                # if there are no galaxies, don't include
                include_cus = False
                
                
            # where to put it?
            if sortIntoFolders:
                if include_cus and include_vir:
                    include_folder = 'associated'
                elif include_cus or include_vir:
                    include_folder = 'ambiguous'
                    print 'include_cus, include_vir = ',include_cus, include_vir
                else:
                    include_folder = 'nonassociated'
            else:
                include_folder = ''
                
            # append this result to the master lists
            masterVirList.append(sorted_virList)
            masterCustomList.append(sorted_cusList)
                
            # include a simple position skyplot?
            if includeSkyPlot:
                # plot
                fig = figure(figsize=(12,10))
                ax = fig.add_subplot(111)
                plot1 = ax.scatter(plotRA,plotDec,marker='d')
                plot2 = ax.scatter(float(AGNposition[0]),float(AGNposition[1]),marker='*')
                ax.set_xlabel('RA')
                ax.set_ylabel('Dec')
        
                savefig('{0}{1}/map2_{2}_{3}_simple.pdf'.format(saveDirectory,include_folder,AGNname,center),format='pdf')

            # save the map plot tables
            if saveMapTables:
                writerOutFile = open('{0}{1}/map_{2}_{3}_table.csv'.format(saveDirectory,include_folder,AGNname,center),'wt')
            
                writer = csv.DictWriter(writerOutFile, fieldnames=fieldnames)
                headers = dict((n,n) for n in fieldnames)
                writer.writerow(headers)
               
                for s in sorted_virList:
                    likelihood, rest = s
                    row = dict((f,o) for f,o in zip(fieldnames,rest))
                    writer.writerow(row)
    
                writerOutFile.close()
            
##########################################################################################
##########################################################################################

            # plot AGN target in center
            ax.scatter(0,0,s=200,c='black',marker='*')
    
#             cbar = plt.colorbar(plot1,ticks=range(0,21),cmap=colmap,orientation='vertical')
#             cbar = plt.colorbar(plot1,ticks=ticks,cmap=colmap,orientation='vertical')
            
            cbar = plt.colorbar(plot1,ticks=ticks,format=r'$\rm %d$',cmap=colmap,orientation='vertical')
            
#             cbar.ax.set_yticklabels([r'$\rm %d$' % i for i ticks])
#             cbar.ax.set_yticklabels(['%d' % i for i ticks])
#             ax.yaxis.set_ticklabels(['%.2f' % 0.1/100*i for i in np.arange(0,100,10)]) 
            cbar.set_label(r'$\rm \Delta v ~[km ~s^{-1}]$')
        
            ax.grid(b=None,which='major',axis='both')
            ax.set_ylim(-500,500)
            ax.set_xlim(-500,500)
            
            ax.set_xlabel(r'$\rm R.A. ~Separation ~[kpc]$')
            ax.set_ylabel(r'$\rm Dec. ~Separation ~[kpc]$')
            if includeTitle:
                title("{0} sightline map velocity = {1} +-/ {2} km/s".format(AGNname,center,velocityWindow))


            # now write it all to file, or display the finished figure
            if saveMaps:
                savefig('{0}{1}/map_{2}_{3}.pdf'.format(saveDirectory,include_folder,AGNname,center),\
                bbox_inches='tight',format='pdf')
            else:
                show()
        
        else:
            print 'TARGET: {0} = NO GALAXIES!!'.format(i)
            
            
#########################################################################################
#########################################################################################
#########################################################################################
#########################################################################################
            
            
    def writeFullResults(listVir,listCustom):
        # write the full results to a new file,  specified by 'outputFile'
        # 
        # this new file will contain all results, both positive and negative
        '''
            list contains entries like [likelihood, [all]], where all is:
            
                ('AGNname',\
                'center',\
                'galaxyName',\
                'environment',\
                'degreesJ2000RA_DecAGN',\
                'degreesJ2000RA_DecGalaxy',\
                'likelihood',\
                'likelihood_1.5'
                'virialRadius',\
                'd^1.5',\
                'impactParameter (kpc)',\
                'redshiftDistances',\
                'vcorrGalaxy (km/s)',\
                'vel_diff',\
                'distGalaxy (Mpc)',\
                'AGN S/N',\
                'majorAxis (kpc)',\
                'minorAxis (kpc)',\
                'inclination (deg)',\
                'positionAngle (deg)',\
                'azimuth (deg)',\
                'RC3flag',\
                'RC3type',\
                'RC3inc (deg)',\
                'RC3pa (deg)',\
                'morphology',\
                'galaxyRedshift')
                
            listVir is sorted by likelihood, listCustom is sorted by custom likelihood
                
        '''
        
        # open the origin file
        f = open(targetFile,'rU')
        reader = csv.DictReader(f)
        
        
        # fieldnames for the new outputFile
        fieldnames = ('AGNname',\
        'center',\
        'galaxyName',\
        'environment',\
        'degreesJ2000RA_DecAGN',\
        'degreesJ2000RA_DecGalaxy',\
        'likelihood',\
        'likelihood_1.5',\
        'virialRadius',\
        'd^1.5',\
        'impactParameter (kpc)',\
        'redshiftDistances',\
        'vcorrGalaxy (km/s)',\
        'radialVelocity (km/s)',\
        'vel_diff',\
        'distGalaxy (Mpc)',\
        'AGN S/N',\
        'majorAxis (kpc)',\
        'minorAxis (kpc)',\
        'inclination (deg)',\
        'positionAngle (deg)',\
        'azimuth (deg)',\
        'RC3flag',\
        'RC3type',\
        'RC3inc (deg)',\
        'RC3pa (deg)',\
        'morphology',\
        'final_morphology',\
        'galaxyRedshift',\
        'AGNredshift',\
        'spectrumStatus',\
        'include',\
        'include_vir',\
        'include_custom',\
        'Lya_v',\
        'vlimits',\
        'Lya_W',\
        'Na',\
        'b',\
        'identified',\
        'comment')
        
        # open the output file
        writerOutFile = open(outputFile,'wt')
        writer = csv.DictWriter(writerOutFile, fieldnames=fieldnames)
        headers = dict((n,n) for n in fieldnames)
        writer.writerow(headers)
        
        count = -1
#         lenList = len(listVir)
        for l in reader:
            AGNname = l['AGNname']
            lyaV = l['Lya_v']

        
            if isNumber(lyaV):
                # this should then correspond to an element in 'list'
                count +=1
                lyaV = float(lyaV)
                
                # grab the corresponding data. list[count] is all the galaxies around a
                # particular target and line
                
                systemsV = listVir[count]
                systemsC = listCustom[count]
                
                # test to make sure this is the correct one. If there were no galaxies,
                # the count will get off
                v1 = systemsV[0][1]
                v1_agnName = v1[0]
                v1_center = v1[1]
                
                if v1_agnName == AGNname and v1_center == lyaV:
                    match = True
                else:
                    match = False
                    count-=1
                
                # now find the 1st, 2nd highest likelihood for both estimates
                sysLen = len(systemsV)
                
                finalEntry = []
                
                if sysLen >=2 and match:
                    # systemsV[0] is the most likely as such: [likelihood, [all]]
                    # thus systemsV[1] is the second most likely
                    vir_1 = systemsV[0]
                    vir_1_like = vir_1[0]
                    vir_1_all = vir_1[1]
                    vir_1_galaxy = vir_1_all[2]
                    
                    vir_2 = systemsV[1]
                    vir_2_like = vir_2[0]
                    vir_2_all = vir_2[1]
                    vir_2_galaxy = vir_2_all[2]
                                    
                    # now for the second, 'custom' likelihood
                    cus_1 = systemsC[0]
                    cus_1_like = cus_1[0]
                    cus_1_all = cus_1[1]
                    cus_1_galaxy = cus_1_all[2]
                    
                    cus_2 = systemsC[1]
                    cus_2_like = cus_2[0]
                    cus_2_all = cus_2[1]
                    cus_2_galaxy = cus_2_all[2]
                    
                    # now decide to marked "include" for each likelihood estimate
                    #
                    # first: for virial radius based likelihood                    
                    if vir_2_like * rigor <= vir_1_like:
                        if vir_1_like >= l_min:
                            virInclude = True
                        else:
                            virInclude = False
                    else:
                        virInclude = False
                    
                    # second: for custom based likelihood                    
                    if cus_2_like * rigor <= cus_1_like:
                        if cus_1_like >= l_min:
                            cusInclude = True
                        else:
                            cusInclude = False                            
                    else:
                        cusInclude = False
                        
                    print
                    print
                    print 'virInclude: ',virInclude
                    print 'vir_1 : ',vir_1
                    print 'vir_2 : ',vir_2
                    
                    print
                    
                    print 'cus_Include: ',cusInclude
                    print 'cus_1: ',cus_1
                    print 'cus_2: ',cus_2
                             
                    
                    # decide on the final thing to include.
                    # first, check if both methods agree that a best galaxy is found
                    if cusInclude and virInclude:
                        
                        # now check if they find the SAME galaxy
                        if vir_1_galaxy == cus_1_galaxy:
                        
                            # in this case vir_1_all and cus_1_all are the same, so just 
                            # pick one to include
                            entry = vir_1_all
                            entry.append(l['AGNredshift'])
                            entry.append(l['spectrumStatus'])
                            entry.append('?')
                            entry.append(virInclude)
                            entry.append(cusInclude)
                            entry.append(l['Lya_v'])
                            entry.append(l['vlimits'])
                            entry.append(l['Lya_W'])
                            entry.append(l['Na'])
                            entry.append(l['b'])
                            entry.append(l['identified'])
                            entry.append(l['comment'])
                            finalEntry = [entry]
                        
                        # if they don't find the same galaxy...
                        else:
                            # enter two lines, first for the virial result, then for 
                            # the custom one later
                            entry = vir_1_all
                            entry.append(l['AGNredshift'])
                            entry.append(l['spectrumStatus'])
                            entry.append('?')
                            
                            # force True for virInclude
                            entry.append(True)
                            # force False for cusInclude
                            entry.append(False)
                            
                            entry.append(l['Lya_v'])
                            entry.append(l['vlimits'])
                            entry.append(l['Lya_W'])
                            entry.append(l['Na'])
                            entry.append(l['b'])
                            entry.append(l['identified'])
                            
                            # update the comment to include this note
                            comment1 = str(l['comment']) + ' : VIRIAL VS CUSTOM RESULT MISMATCH'
                            entry.append(comment1)
                            
                            
                            # now enter a second line for the custom result
                            entry2 = cus_1_all
                            entry2.append(l['AGNredshift'])
                            entry2.append(l['spectrumStatus'])
                            entry2.append('?')
                            
                            # now force False for virInclude
                            entry2.append(False)
                            # and True for cusInclude
                            entry2.append(True)
                            
                            entry2.append(l['Lya_v'])
                            entry2.append(l['vlimits'])
                            entry2.append(l['Lya_W'])
                            entry2.append(l['Na'])
                            entry2.append(l['b'])
                            entry2.append(l['identified'])
                            
                            # include the same updated comment as above
                            entry2.append(comment1)
                            
                            finalEntry = [entry,entry2]
                            
                    # if the methods do not agree (i.e. one finds a best galaxy, and the
                    # other finds NO best galaxy)
                    else:
                        # in this case vir_1_all and cus_1_all will not be the same, so I
                        # need to pick the right one
                        if virInclude:
                            entry = vir_1_all
                        else:
                            entry = cus_1_all
                            
                        entry.append(l['AGNredshift'])
                        entry.append(l['spectrumStatus'])
                        entry.append('?')
                        entry.append(virInclude)
                        entry.append(cusInclude)
                        entry.append(l['Lya_v'])
                        entry.append(l['vlimits'])
                        entry.append(l['Lya_W'])
                        entry.append(l['Na'])
                        entry.append(l['b'])
                        entry.append(l['identified'])
                        
                        # update the comment to say the includes do not match
                        comment1 = str(l['comment']) + " : VIRIAL AND CUSTOM INCLUDE MISMATCH"
                        entry.append(comment1)
                        
                        finalEntry = [entry]              
                            
                        
                
                # now if there is only one galaxy
                elif sysLen == 1 and match:
        
                    # is the likelihood big enough? (>= l_min)
                    vir_1 = systemsV[0]
                    vir_1_like = vir_1[0]
                    if float(vir_1_like) >= l_min:
                        virInclude = True
                    else:
                        virInclude = False
                
                    cus_1 = systemsC[0]
                    cus_1_like = cus_1[0]
                    if float(cus_1_like) >= l_min:
                        cusInclude = True
                    else:
                        cusInclude = False
                    
                    entry = vir_1[1]
                    entry.append(l['AGNredshift'])
                    entry.append(l['spectrumStatus'])
                    entry.append('?')
                    entry.append(virInclude)
                    entry.append(cusInclude)
                    entry.append(l['Lya_v'])
                    entry.append(l['vlimits'])
                    entry.append(l['Lya_W'])
                    entry.append(l['Na'])
                    entry.append(l['b'])
                    entry.append(l['identified'])
                
                    # include the same updated comment as above
                    comment1 = str(l['comment']) + " : ONLY 1 GALAXY"
                    entry.append(comment1)
                
                    finalEntry = [entry]

                # no galaxies, mark it so in the comments
                else:
                    virInclude = False
                    cusInclude = False
                    
                    # also finalInclude is forced to FALSE - see below

                    comment1 = str(l['comment']) + " : NO GALAXIES"
                    
                    finalEntry = [[l['AGNname'],\
                    lyaV,\
                    'x',\
                    0,\
                    l['degreesJ2000RA_DecAGN'],\
                    ('x','x'),\
                    -99,\
                    -99,\
                    -99,\
                    -99,\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    l['AGN S/N'],\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    'x',\
                    l['AGNredshift'],\
                    l['spectrumStatus'],\
                    False,\
                    virInclude,\
                    virInclude,\
                    l['Lya_v'],\
                    l['vlimits'],\
                    l['Lya_W'],\
                    l['Na'],\
                    l['b'],\
                    l['identified'],\
                    comment1]]
                    
                                        
                # write it all out to file
                for entry in finalEntry:
#                     print 'entry: ',entry

                    row = dict((f,o) for f,o in zip(fieldnames,entry))
                    writer.writerow(row)
            
            
        # close the files
        f.close()
        writerOutFile.close()
            
            
##########################################################################################
##########################################################################################
       
    # now actually do it
    if saveResults:
        writeFullResults(masterVirList,masterCustomList)
    
    print "Finished."
    

if __name__=="__main__":
    main()