# Open a tvstudy.txt file, find stations that fail the interference check, and 
# determine power reduction required to bring interference below de minimis amount.
# Generate a kml file showing interference
# Load pattern from TVStudy and show modified pattern after interference analysis
# Created by Doug Lung - dlung@transmitter.com
# Version 1.2 - provide pattern plot in text file and TVStudy format
# Version 1.3 - get call letters from sources.csv and add call letters to kml ploy
# Version 1.4 - Fix crash when there is no comment and when **IX check failure appears before Scenario comparisons
# Version 1.5 - Fix error when non-directional antenna pattern is used
# Version 1.6 - Generate patterns with 36 and 72 degree radials - plot and tabulate in csv file for TVStudy import
# Version 1.7 - Add support for rotated antenna patterns and odd azimuths exceeding the number allowed on old FCC forms
# Version 1.8 - Find correct pattern data if more than desired station in parameters.csv file
# Version a2.0 - Alpha version - first attempt to modify program to work with TVStudy 2.1 - not functional
# Version 2.0a1 - Appears to be working!!
# Version 2.0a2 - Modified parameters.csv file import to match TVStudy 2.2 format
# Version 2.0a3 - Fix bug with LPTV protection of full power stations to use 0.5% threshold
# Version 2.0a3 - Create sources list from individual study directories instead of "coverage' or 'MX#1' 
# Version 2.0a3 - Appears to work okay with TVStudy 2.2.2 May update
# Version 2.0a4 - Add DD (DTS) stations in those to study with 0.5% interference threshold
# Version 2.0a4 - Fix parsing of parameters.csv antenna pattern to ignore \n and space characters
# Version 2.0a5 - Use sources.csv in Coverage directory to determine SourceID - fix issue with replicated stations
# Version 2.0a6 - Correct pattern extraction offset error, use first SourceID if multiple ones in sources.csv
# Version 2.0a7 - Change matplotlib depreciated "bg" to "facecolor" to eliminate error message. 
#    Added test for mutually exclusive applications **MX as well as **IX
# Version 2.0a8 - Added option with StudyScenario and IgnoreScenario to only test certain scenarios. 
#    Added line with call, filenumber and scenario number to standard output amd separated scenarios with blank line
# Version 2.0a9 - Fix crash in pattern plotting functions when using updated matplotlib
# Version 2.0a10 - Modified to work with tvixstudy.txt changes in TVStudy 2.2.4 (Use 2.0a9 for TVStudy 2,2,3) 
# Version 2.0a11 - Modifed code to work with tvixstudy.txt from either TVStudy 2.2.4 or TVStudy 2.2.3
# Version 2.0a12 - Use study file number instead of Fac ID to find interference source in sources.csv (2.2.4 changed order in sources.csv)
#    2.0a12 should work with both TVStudy 2.2.3 and 2.2.4 output files
# Version 2.0a13 - Fix variable conflict when "StudyScenario" and "IgnoreScenario" are used.
# Version 2.0a14 - Modify tvixstudy.txt parsing to ignore Proposal "before" file number when scenarios are edited
# Version 2.0a15 - Truncate long file numbers at 22 characters when doing sources.csv lookup
# Version 2.0a16 - Modify function for parsing parameters.csv file to eliminate crash caused when city name contains a comma
# Version 2.1a0 - TVStudy 2.3 fixes
# Version 2.1a2 and 2.1a3 - More TVStudy 2.3 fixes with tvixstudy.txt parsing. Changed pattern plots to png to avoid matplotlib error
# Version 2.1a4 - Modify routines to only consider interference in a desired station's own country 
# Version 2.1a5 - Modified test at Line 165 - 167 to allow for added callsign in tvstudy output
# Version 2.1a6 - Change "Proposal" directory to "Proposal_NoIX"
# Version 2.1b1 - Incorporate geo.py code in program
# Version 2.1b2 - Test for simplekml and matplotlib modules. If not available, warn and disable functions requiring them
# Version 2.1b3 - Provide text output in 'ixstudy.txt' file
# Version 2.1b4 - Modify to handle directory name changes in TVStudy 2.3.1 and cover MX cases 

import os
import math
from operator import itemgetter
try:
    import simplekml
    UseSimpleKML = True
except: 
    UseSimpleKML = False
    print("")
    print("simplekml module is not installed or out-of-date. Map creation disabled.")
    print("")
# import geo
import numpy as np
try:
    import matplotlib.pyplot as plt
    import matplotlib as mpl
    UseMatplotlib = True
except:
    UseMatplotlib = False
    print("")
    print("matplotlib not installedd or out-of-date. Pattern graphs disabled.")
    print("")              

# Initial routine to test only specific scenarios
# Scenarios are identified by file number and scenario number
# When a entries are present, the interference test routine will see if the filenumber is in list
# and if so, if the scenario number is in the list before adding the study to the 
# ixList for StudyScenario or removing it from the list for IgnoreScenario. 
# If the scenario number is "ALL", all scenarios with the filenumber will be studied or
# ignored, depending on whether it is in the StudyScenario or IgnoreScenario dictionary. . 
StudyScenario = {}  # Do not delete this line - add entries below
# StudyScenario["ixFileNum"] = ['1','2','3']
# Examples - use "ALL" to study all scenarios for the filenumber. 
# Note all items are strings and must have quotes, even scenario numbers
# StudyScenario["DTVBL13602"] = ['1','2','3','4']
# StudyScenario["DTVBL13602"] = "ALL"

IgnoreScenario = {} # Do not delete this line - add entries below
# IgnoreScenario[ixFileNum] = ['1','2','3']
# Examples - use "ALL" to ignore all scenarios for the filenumber. 
# Note all items are strings and must have quotes, even scenario numbers
# IgnoreScenario["DTVBL13602"] = ['2','3','4']
# IgnoreScenario["DTVBL13602"] = "ALL"
IgnoreScenario["BLANK0000293163"] = "ALL"

# The functions below are from geo.py and provide angle calculations
# geo.py is Copyright (C) 2010  Maximilian Hoegner <hp.maxi@hoegners.de>
# Distributed under terms of GNU General Public License Version 3 or later
# See <http://www.gnu.org/licenses/>.


EARTH_RADIUS = 6370000.
MAG_LAT=82.7
MAG_LON=-114.4

def xyz(lat,lon,r=EARTH_RADIUS):
    """ Takes spherical coordinates and returns a triple of cartesian coordinates """
    x = r*math.cos(math.radians(lat))*math.cos(math.radians(lon))
    y = r*math.cos(math.radians(lat))*math.sin(math.radians(lon))
    z = r*math.sin(math.radians(lat))
    return x,y,z

def cross(p1,p2):
    """ Cross product of two vectors """
    x = p1[1]*p2[2]-p1[2]*p2[1]
    y = p1[2]*p2[0]-p1[0]*p2[2]
    z = p1[0]*p2[1]-p1[1]*p2[0]
    return x,y,z

def determinant(p1,p2,p3):
    """ Determinant of three vectors """
    return dot(p1,cross(p2,p3))

def dot(p1,p2):
    """ Dot product of two vectors """
    return p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2]

def normalize_angle(angle):
    """ Takes angle in degrees and returns angle from 0 to 360 degrees """
    cycles = angle/360.
    normalized_cycles = cycles - math.floor(cycles)
    return normalized_cycles*360.

def sgn(x):
    """ Returns sign of number """
    if x==0: return 0.
    elif x>0: return 1.
    else: return -1.

def angle(v1,v2,n=None):
    """ Returns angle between v1 and v2 in degrees. n can be a vector that points to an observer who is looking at the plane containing v1 and v2. This way, you can get well-defined signs. """
    if n==None:
        n=cross(v1,v2)
    prod = dot(v1,v2) / math.sqrt( dot(v1,v1) * dot(v2,v2) )
    if prod>1: prod=1.0 # avoid numerical problems for very small angles
    rad = sgn(determinant(v1,v2,n)) * math.acos( prod )
    deg = math.degrees(rad)
    return normalize_angle(deg)

def great_circle_angle(p1,p2,p3):
    """ Returns angle w(p1,p2,p3) in degrees. Needs p1 != p2 and p2 != p3. """
    n1=cross(p1,p2)
    n2=cross(p3,p2)
    return angle(n1,n2,p2)

geographic_northpole=xyz(90,0)
# End of functions from geo.py


# Create and open "ixstudy.txt" file for study output
ixstudyFile = 'ixstudy.txt' 
ixstudy = open(ixstudyFile,'w')
global ixstudyText
ixstudyText = ""

# This function looks for a tvixstudy.txt file in the current directory 
# It takes no arguments
# Returns two lists:
#   StudyInfo = Study Name, Comment, CallSign, TxLatitude, TxLongitude, PeakERP, SourceFacID
#   IxList, for each station: ixCall, ixStatus, ixScenarioNum, ExistingIxFreePop, StudyIxFreePop, i'xService, ixFileNumber
# Tested and working with TVStudy 2.21

# No complete sources.csv exists in Version 2.2.2 - create one and update as U_CelData.csv files are loaded

sourcesDict = {}

# In most cases the source ID to study will be "1", but if a station is replicated it could be "2"
# This function opens the Coverage/sources.csv file for the studied facility and returns all source IDs
# in this file in a list. For now, only the first line is used as the source ID to study.
def getSourceID(SourceFileNum):
    global ixstudyText
    try:
        sourceIDcsv = open('Proposal_NoIX/cellcsv/sources.csv','r')
    except:
        print("Proposal_NoIX/cellcsv/sources.csv file not found. Exiting")
        exit()
    sourceIDlist = []
    if len(SourceFileNum) > 22 :
        SourceFileNum = SourceFileNum[0:22]
    print('Looking for File number ',SourceFileNum)
    for line in sourceIDcsv:
        print(line)
        ixstudyText = ixstudyText + (line) + "\n"
        if (line.strip().find(SourceFileNum) > 0):
            sourceIDlist.append(line.split(',')[0])
    if(sourceIDlist == []) :
        print("SourceID for File Number ", SourceFileNum, " not found. Exiting")
        exit()
    return(sourceIDlist)
                
# This function parses the tvixstudy.txt file to get information on the station being studied and
# stations receiving interference over de minimis levels.
def getTVstudyData() :
    global ixstudyText
    try:
        tvixstudy = open('tvixstudy.txt','r')
    except:
        print("tvixstudy.txt file not found. Exiting")
        exit()
    
    Comment = 'No comment'
    SourceFileNumber = ''
    SourceFacID = ''
    ScenarioComparisons = False
    Interference_test = False
    Desired_Found = False
    Desired_Data = False
    NextLine = ''
    ixList = []
    for line in tvixstudy:
# Change find index comparison for TVStudy 2.3.0
        if (line.strip().find('Study:') != -1):
            lineList = line.split(':')
# Change index to "1" for TVStudy 2.3.0
            StudyAll = lineList[1].strip()
            Study = StudyAll.split(',')[0]
        if (NextLine == 'Comment'):
            Comment = line.strip()
            NextLine = ''
        if (line.strip().find('Comment:') >= 0):
            NextLine = 'Comment'
        if (line.strip().find('Proposal:') >= 0):
            lineList = line.strip().split(':')
            CallSign = lineList[1].strip().split()[0]
        if (line.strip().find('Facility ID:') >= 0) :
            lineList = line.strip().split(':')
            if (SourceFacID == '') : # Don't change source facility ID if has already been stored
                SourceFacID = lineList[1].strip().split()[0]
        if (line.strip().find('File number:') >= 0) :
            lineList = line.strip().split(':')
            if (SourceFileNumber == '') :  # Don't change file number if it has already been stored
                SourceFileNumber = lineList[1].strip().split()[0]
                print('SourceFileNumber = ',SourceFileNumber)
                ixstudyText = 'SourceFileNumber = ' + SourceFileNumber
        if (line.strip().find('Latitude: ') >= 0):
            lineList = line.strip().split(':')
            LatDeg = int(lineList[1].split()[0])
            LatMin = int(lineList[1].split()[1])
            LatSec = float(lineList[1].split()[2])
        if (line.strip().find('Longitude: ') >= 0):
            lineList = line.strip().split(':')
            LonDeg = int(lineList[1].split()[0])
            LonMin = int(lineList[1].split()[1])
            LonSec = float(lineList[1].split()[2])
        if (line.strip().find('Peak ERP:') >= 0):
            lineList = line.strip().split(':')
            PeakERP = lineList[1].strip().split()[0]
        if ((line.strip().find('**IX:') >= 0) | (line.strip().find('**MX:') >= 0)) & Interference_test :
            ScenarioComparisons = True
            Interference_test = False
        if (line.strip().find('Interference to ') >= 0) :
            lineList = line.strip().split('scenario')
# Modified to account for added callsigns in TVStudy 2.3
            if len(lineList[0].strip().split()) == 5 :
                Interference_test = True
                ixFileNumber = lineList[0].strip().split()[4]
                ixStatus = lineList[0].strip().split()[3]
                if(ixStatus.strip().find(',') >= 0) :
                    ixStatus = ixStatus.strip().split(',')[0]
                ixScenarioNum = lineList[1].strip().split()[0]
            else :
                Interference_test = False
        if (line.strip().find('Desired') >= 0) & ScenarioComparisons :
            lineList = line.strip().split()
            ixCall = line[13:23].strip()
            ixService = line[29:31]
#            print(ixCall + ' ' + ixScenarioNum)
            ScenarioComparisons = False
            Desired_Found = True
        if Desired_Data :
            lineList = line.strip().split()
            ExistingIxFreePop = lineList[5].replace(',','')
            StudyIxFreePop = lineList[7].replace(',','')
            # Determine country for service and interference populations)
            if len(lineList) > 10 :
                CountryList = line.strip().split('(')
                ixCountry = CountryList[1].strip().split('in ')[1]
                ixCountry = ixCountry[0:-1]
            else : 
                ixCountry = ''
                
            if ixCall[0:1] == 'X' :
                srcCountry = '3'
            elif ixCall[0:1] == 'C' :
                srcCountry = '2'
            else : 
                srcCountry = '1'
            
            # IxCall, status, scenarioNum, ExistingIxFreePop, StudyIxFreePop, ixService, ixFileNumber
            if ixCountry == '' :
                ixList.append([ixCall, ixStatus, ixScenarioNum, ExistingIxFreePop, StudyIxFreePop, ixService, ixFileNumber, srcCountry])
                Desired_Data = False
        if (line.strip().find('Service area') >= 0) & Desired_Found :
            Desired_Data = True
            Desired_Found = False
    
    TxLatitude = float(LatDeg) + float(LatMin)/60 + float(LatSec)/3600
    TxLongitude = float(LonDeg) + float(LonMin)/60 + float(LonSec)/3600
    PeakERP = float(PeakERP)
    StudyInfo = [Study, Comment, CallSign, TxLatitude, TxLongitude, PeakERP, SourceFacID, SourceFileNumber]
    tvixstudy.close()
    
    for i, ixCase in enumerate(ixList):
        ixFileNumber = ixCase[6]
        ixStatus = ixCase[1]
        ixScenario = ixCase[2]
        #Modified for TVStudy 2.3 and 2.3.1
# Remove "#" from scenario number to work with TVStudy 2.3.1 directories
# Don't try to find a U_CelData.csv file in the scenario for the facility being studied
        if SourceFileNumber != ixFileNumber:
            U_CelDataFile = "./IX_" + ixFileNumber + "_" + ixStatus + "_" + ixScenario + "_after/cellcsv/U_CelData.csv"
            try :
                U_CelData = open(U_CelDataFile, 'r')
            except :
                print('Unable to open ', U_CelDataFile, ' exiting')
                exit()
            U_CelData.close()
    return(StudyInfo, ixList)
  
# This function loads the TVStudy parameters.csv file for the station under study
# and extracts the antenna pattern.
# The antenna pattern is interpolated to 361 radials 
# The function returns a dictionary with the 361 degree patterns
# The "0" radial has to be repeated in the dictonary as a "360' degree azimuth to allow
# and graphing to work properly. 
# Tested and working with TVstudy 2.2
# If erroneous pattern data is retreived, check the offsets used in creating and extracting
# data from the PatternLine string.
def getTVstudyPattern(SourceFacID) :
    try:
        #Modified for TVStudy 2.3
        patternCsv = open('Proposal_NoIX/parameters.csv','r')
    except:
        print("parameters.csv file not found. Exiting")
        exit()
    DataLine = 0
    patternLine = []
    finalPatternLine = []
    for item in patternCsv:
        if DataLine & (len(patternLine) < 1) :
            patternLine = item.strip().split(',')
            if (patternLine[0].strip()) != SourceFacID :
                patternLine = []
        if ('FacID' in item) :
            DataLine = 1
# Look for quoted strings with commas in them and replace comma with space
    count = 0
    for item in patternLine :
        if len(item) > 0 :  # Test for zero length items
            if ((item[0] != '"') & (item[-1] != '"')) | ((item[0] == '"') & (item[-1] == '"')): # Don't mess with items without double quotes
                finalPatternLine.append(item)
            elif (item[0] == '"') & (item[-1] != '"') :  # If an open double quote is found without a closing quote, temporarily save it
                tempA = item
                tempB = ''
            elif (item[0] != '"') & (item[-1] == '"') : # When the closing quote is found, combine with the segment saved above and add to list
                tempB = item
                finalPatternLine.append(tempA + ' ' + tempB)
                tempA = ''
                tempB = ''
        else :
            finalPatternLine.append(item) # Save the zero length items to preserve the index values and avoid dropping azimuths
         
    patternLine = finalPatternLine # Overwrite patternLine list with finalPatternLine list to avoid rewriting the rest of the function
    pLlength = len(patternLine)
#    print('PATTERN LINE: ',patternLine)
    patrotstr = patternLine[22].strip()
#    print(pLlength, patrotstr)
    if patrotstr != '':
        AntRotation = int(round(float(patrotstr),0))
    else :
        AntRotation = 0
#    print(pLlength, patrotstr)
    if len(patternLine) > 40 :
        pattern10 = patternLine[40:76]
#        print(pattern10)
        patternDict = {}
        angle = 0 + AntRotation
        for item in pattern10 :
#            if item == '' :
#                rf = 1.0
#            if (item != '\n') & (item != ' ') :
            if item != '' :
#                print(item)
                rf = float(item)
            patternDict[angle] = rf
            angle = angle + 10
            if angle >= 360 :
                angle = angle - 360

        patternOdd = patternLine[77:]
        angle = ''
        for item in patternOdd :
            if item != '':
                if angle == '' : 
                    angle = int(item) + AntRotation
                    if angle >= 360 : 
                        angle = angle - 360
                else:
                    rf = float(item.strip())
                    patternDict[angle] = rf
                    angle = ''

        AzSorted = sorted(patternDict.keys())
        RfSorted = []
        for item in AzSorted :
            RfSorted.append(patternDict[item])
        RfSorted.append(patternDict[AzSorted[0]])
        AzSorted.append(360)
    else :
        AzSorted = [0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,\
            210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360]
        RfSorted = [1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,\
            1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]

    Az360 = range(361)
    Rf360 = []
    Pattern360 = {}
    for Azimuth in Az360 :
        Rf = (np.interp(Azimuth, AzSorted, RfSorted))
        Rf360.append(Rf)
        Pattern360[Azimuth] = Rf
    return(Pattern360)
            
# This function does a cell by cell interference analysis 
# Arguments:
#   SourceID = Source ID of station being studied - undesired in this case
#   StudyInfo, IxData (list from IxList)
#   StudyInfo = [Study Name, Comment, CallSign, TxLatitude, TxLongitude, PeakERP, SourceFacID]
#   IxData = [ixCall, ixStatus, ixScenarioNum, ExistingIxFreePop, StudyIxFreePop, ixService, ixFileNumber,srcCountry] 
# Removed threshold as it can be calculated from ixService in IxData
# Returns:
#   List of study results: VictimCall, scenario, ExistingIxFreePop, StudyIxFreePop, IxPopOverLimit, IxPercent, ERPreductionReq, ixThreshold
#   List of cells to clear: CellKey, PointLat, PointLon, PointPop, IxLevelDB, BearingToTx
# Generates: CSV files for ix from source and cells to clear, kml files for unique ix from source and cells to clear 
def VictimIxStudy(SourceID, StudyInfo, IxData) :
    global ixstudyText
    ixCall = IxData[0]
    ix_source = SourceID
    ixStatus = IxData[1]
    ixScenario = IxData[2]
    ixService = IxData[5]
    ixFileNumber = IxData[6]
    srcCountry = IxData[7]
    if (ixService == 'DT') | (ixService == 'DC') | (ixService == 'DD') :
        ixThreshold = 0.005
    else :
        ixThreshold = 0.02
# Remove "#" from scenario number to work with TVStudy 2.3.1    
    CsvFile = "./IX_" + ixFileNumber + "_" + ixStatus + "_" + ixScenario + "_after/cellcsv/U_CelData.csv"
# Remove "#" from scenario number to work with TVStudy 2.3.1
    SourcesFile = "./IX_" + ixFileNumber + "_" + ixStatus + "_" + ixScenario + "_after/cellcsv/sources.csv"
    try:
        sourcesCsv = open(SourcesFile,'r')
    except:
        print("sources.csv file not found. Exiting")
        exit()
    for item in sourcesCsv :
        source = item.split(',')
        key = source[0].strip()
        facID = source[8].strip()
        fileNum = source[10].strip()
        call = source[11].strip()
#        print("Source: ",facID + ' ' + fileNum + ' ' + call)
        sourcesDict[key] = [facID, fileNum, call]
    try:
        csvData = open(CsvFile, mode='r')
        print('Opened ',CsvFile)
    except:
        print(CsvFile + ' not found. Exiting')
        exit()
    Current_ix_free_pop = int(IxData[3])
    Study_ix_free_pop = int(IxData[4])
    study_ERP = StudyInfo[5]
    sourceCall = StudyInfo[2]
    count = 0
    VictimIxData = []
    data = []

    # Save only cells that have both interference and population
    # Added comparison with srcCountry to count only interference in station's country
    for line in csvData :
        if (line.strip() != ''):
            data = line.strip().split(',')
            cellCountry = data[0].split('-')[2]
        # Only count interference in source country
            if((data[18].strip() == '1') & (int(data[4]) > 0)) & (cellCountry == srcCountry):
                VictimIxData.append(data)
                count = count + 1
    print(ixCall + ' File Number: ' + ixFileNumber + '    Scenario: ' + ixScenario)
    ixstudyText = ixstudyText + ixCall + ' File Number: ' + ixFileNumber + '    Scenario: ' + ixScenario + "\n"
    print("Service = ",ixService)  # DEBUG LINE
    ixstudyText = ixstudyText + "Service = " + ixService + "\n"
    print('Number of cells with interference and population: ', count)
    ixstudyText = ixstudyText + 'Number of cells with interference and population: ' + str(count) + "\n" 
    print('ix_source = ',ix_source)
    ixstudyText = ixstudyText + 'ix_source = ' + ix_source + "\n"
    
    # Determine cells with unique interference from ix_source and count total amount of unique interference
    # You can check this number against TVStudy output to verify correct processing

    VictimIxDataSorted = sorted(VictimIxData,key=itemgetter(0))
    
    ixPop = 0
    data = []
    tempItem = []
    lastCellID = ''
    IxFromSource = []
    for item in VictimIxDataSorted :
        if item[0] != lastCellID :
            if tempItem != [] :
                IxFromSource.append(tempItem)
                ixPop = ixPop + int(tempItem[4])
                tempItem = []
            if item[6] == ix_source : 
                tempItem = item
                tempItem.append((float(tempItem[16]) + float(tempItem[17]) - float(tempItem[15])))
            lastCellID = item[0]
        else :
            tempItem = []
            lastCellID = item[0]

    if(len(tempItem) > 0) :
        IxFromSource.append(tempItem)
        ixPop = ixPop + int(tempItem[4])

    
    print('Number of cells with interference from source: ',len(IxFromSource))
    ixstudyText = ixstudyText + 'Number of cells with interference from source: '+ str(len(IxFromSource)) + '\n'

    print('ixPop from source = ',ixPop)
    ixstudyText = ixstudyText + 'ixPop from source = ' + str(ixPop) + '\n'

    #Sort IxFromSource by calculated, appended IX_level (low to high)
    IxFromSourceSorted = sorted(IxFromSource,key=itemgetter(21))

    # Output sorted data as a CSV file
    ixOutputFile ="Ix_from_" + sourceCall + "_to_" + ixCall + "_" + ixStatus + "_" + ixScenario+'.csv'
    ixOutput = open(ixOutputFile,'w')
    header = 'Key'+','+'Point_Lat'+','+'Point_Lon'+','+'Point_Area'+','+'Point_Pop'+','+'Point_HH'+','+'U_Source_Key'+','+'U_Site_Num'+','+'D_Source_Key'+','+'D_Site_Num'+','+'D_Field_Str'+','+'Serv_Flag'+','+'U_Field_Str'+','+'Rcv_Pat_Adj'+','+'R2T_Bearing'+','+'DU_Ratio'+','+'DU_Thres'+','+'SNR_Fact'+','+'Causes_Ix'+','+'D_KWX_Flag'+','+'U_KWX_Flag' + ',' + 'IX_level \n'
    ixOutput.write(header)
    for line in IxFromSourceSorted :
        csvout = line[0] + ',' + line[1] + ',' + line[2] + ',' + line[3] + ',' + line[4] + ',' + line[5] + ',' + line[6] + ',' + line[7] + ',' + line[8] + ',' + line[9] + ',' + line[10] + ',' + line[11] + ',' + line[12] + ',' +  line[13] + ',' + line[14] + ',' + line[15] + ',' + line[16] + ',' + line[17] + ',' + line[18] +  ',' + line[19] +  ',' + line[20] +  ',' +  '{:2.3f}'.format(float(line[21])) + '\n'      
        ixOutput.write(csvout)

    ixOutput.close
    csvout = ''

    # Calculate de minimis interference and round up to next integer. Add to existing unique interference to determine new interference limit
    ixNewAllowed = math.floor(Current_ix_free_pop * ixThreshold) # math.floor rounds down result to nearest integer
    ixOverLimit = Current_ix_free_pop - Study_ix_free_pop - ixNewAllowed

    print('Current_ix_free_pop = ' + str(Current_ix_free_pop))
    ixstudyText = ixstudyText + 'Current_ix_free_pop = ' + str(Current_ix_free_pop) + "\n"
    print('de minimis threshold = ' + '{:2.3f}'.format(float(ixThreshold)))
    ixstudyText = ixstudyText + 'de minimis threshold = ' + '{:2.3f}'.format(float(ixThreshold)) + "\n"
    print('ixNewAllowed = ' + str(ixNewAllowed))
    ixstudyText = ixstudyText + 'ixNewAllowed = ' + str(ixNewAllowed) + "\n"
    print('ixOverLimit = ' + str(ixOverLimit))
    ixstudyText = ixstudyText + 'ixOverLimit = ' + str(ixOverLimit) + "\n"

    # Step through IxFromSourceSorted until interference drops below de minimis IX_level
    popToClear = ixOverLimit
    popRemoved = 0
    powerReduction = 0
    CellsToClear = []
    CellsIxPermitted = []

    for cell in IxFromSourceSorted :
        if popRemoved < popToClear :
            Pop = int(cell[4])
            excessPower = cell[21]
            CellsToClear.append(cell)
            popRemoved = popRemoved + Pop
            if excessPower > powerReduction :
                powerReduction = excessPower
        else :
            CellsIxPermitted.append(cell)

    print('Power reduction required = ' + '{:2.3f}'.format(powerReduction) + ' dB')
    ixstudyText = ixstudyText + 'Power reduction required = ' + '{:2.3f}'.format(powerReduction) + ' dB' + "\n"
    maximumERP = study_ERP * (math.pow(10,-(powerReduction/10)))
    print('Maximum ERP allowed for this scenario = ' + '{:4.3f}'.format(maximumERP))
    ixstudyText = ixstudyText + 'Maximum ERP allowed for this scenario = ' + '{:4.3f}'.format(maximumERP) + "\n\n"
    print('\n')

    # Create new csv file with only cells to cells_to_clear
    ixOutputFile =sourceCall + "_ix_cells_to_clear_to_protect_" + ixCall + "_" + ixStatus + "_" + ixScenario+'.csv' 
    ixOutput = open(ixOutputFile,'w')
    header = 'Key'+','+'Point_Lat'+','+'Point_Lon'+','+'Point_Area'+','+'Point_Pop'+','+'Point_HH'+','+'U_Source_Key'+','+'U_Site_Num'+','+'D_Source_Key'+','+'D_Site_Num'+','+'D_Field_Str'+','+'Serv_Flag'+','+'U_Field_Str'+','+'Rcv_Pat_Adj'+','+'R2T_Bearing'+','+'DU_Ratio'+','+'DU_Thres'+','+'SNR_Fact'+','+'Causes_Ix'+','+'D_KWX_Flag'+','+'U_KWX_Flag' + ',' + 'IX_level \n'
    ixOutput.write(header)
    for line in CellsToClear :
        csvout = line[0] + ',' + line[1] + ',' + line[2] + ',' + line[3] + ',' + line[4] + ',' + line[5] + ',' + line[6] + ',' + line[7] + ',' + line[8] + ',' + line[9] + ',' + line[10] + ',' + line[11] + ',' + line[12] + ',' +  line[13] + ',' + line[14] + ',' + line[15] + ',' + line[16] + ',' + line[17] + ',' + line[18] +  ',' + line[19] +  ',' + line[20] +  ',' +  '{:2.3f}'.format(float(line[21])) + '\n'
        ixOutput.write(csvout)
    ixOutput.close
    cvsout = ''

    # Create new kml file showing ix cells and cells to clear
    if (UseSimpleKML == True) :
        kml = simplekml.Kml(open=1)
        for Key, Point_Lat, Point_Lon, Point_Area, Point_Pop, Point_HH, U_Source_Key, U_Site_Num, D_Source_Key, D_Site_Num, D_Field_Str, Serv_Flag, U_Field_Str, Rcv_Pat_Adj, R2T_Bearing, DU_Ratio, DU_Thres, SNR_Fact, Causes_Ix, D_KWX_Flag, U_KWX_Flag, IX_level in CellsToClear :
            pnt = kml.newpoint()
            pnt.name = Point_Pop
            pnt.iconstyle.color = simplekml.Color.red
            pnt.description = 'Point_Pop: ' + Point_Pop + '   DU_Ratio: ' + DU_Ratio + '   DU_Thres: ' + DU_Thres + '   IX_level(dB): ' + '{:2.3f}'.format(float(IX_level)) 
            pnt.coords = [('-'+Point_Lon, Point_Lat)]
        for Key, Point_Lat, Point_Lon, Point_Area, Point_Pop, Point_HH, U_Source_Key, U_Site_Num, D_Source_Key, D_Site_Num, D_Field_Str, Serv_Flag, U_Field_Str, Rcv_Pat_Adj, R2T_Bearing, DU_Ratio, DU_Thres, SNR_Fact, Causes_Ix, D_KWX_Flag, U_KWX_Flag, IX_level in CellsIxPermitted :
            pnt = kml.newpoint()
            pnt.name = Point_Pop
            pnt.iconstyle.color = simplekml.Color.lightgreen
            pnt.description = 'Point_Pop: ' + Point_Pop + '   DU_Ratio: ' + DU_Ratio + '   DU_Thres: ' + DU_Thres + '   IX_level(dB): ' + '{:2.3f}'.format(float(IX_level)) 
            pnt.coords = [('-'+Point_Lon, Point_Lat)]
        kml.save(sourceCall +"_cells_with_ix_to_" + ixCall + "_" + ixStatus + "_" + ixScenario+ ".kml")
 
    sourcesCsv.close()
    csvData.close()
    return(ixCall, Current_ix_free_pop, Study_ix_free_pop, powerReduction, CellsToClear)
    
# Combines results from VictimIxStudy results to generate kml file showing all cells to clear
# Combines required power reductions and bearings to clear all ix cases
def CombineVictims(IxStudies, txLat, txLong) :
    Radials = {}
    ReductionList = []
    CombinedCellsToClear = []
    
#   Iterate list - calculate bearing from transmitter to cell
#   Create dictionary with bearing rounded to nearest degree as key and power reduction as azimuth
#   if same key, save largest power reduction required
    for Cell in IxStudies :
        cellLat = float(Cell[1])
        cellLong = -(float(Cell[2]))                     
        CellLoc = xyz(cellLat, cellLong)
        TxLoc = xyz(txLat,-(txLong))
        bearing = great_circle_angle(CellLoc,TxLoc,geographic_northpole)
        dBreduction = Cell[21]
        ReductionList.append(dBreduction)
        bearingKey = int(round(bearing,0)) # Round to nearest degree

        if bearingKey in Radials :
            if dBreduction > Radials[bearingKey] :
                Radials[bearingKey] = dBreduction
                CombinedCellsToClear.append(Cell)
        else :
            Radials[bearingKey] = dBreduction
            CombinedCellsToClear.append(Cell)
    if ReductionList != [] :
        maxReduction = max(ReductionList)
    else :
        maxReduction = 0
    return(Radials, CombinedCellsToClear, maxReduction)

# Function to display antenna azimuth patterns on a polar plot
def plotAzimuthPattern(PatternDict1, plotFile) :
    angles = np.arange(0,361.0,1)
    Pattern1 = []
    for radial in angles :        # Put relative fields in a list start with 0 degree radial
        Pattern1.append(PatternDict1[radial])        
    theta = (np.pi/180.0) * angles
    mpl.rc("ytick",labelsize = "small")
    fig1 = plt.figure(figsize = (5,5.4), dpi=120, facecolor='White')
    ax1 = fig1.add_axes([0.1,0.1,0.8,0.8], polar=True)
    ax1.set_theta_zero_location("N")
    ax1.set_theta_direction(-1)
    ax1.set_ylim(0,1)
    ax1.set_yticks(np.arange(0,1,0.1))
    ax1.plot(theta,Pattern1,lw=2.5)
    ax1.set_title(plotFile, position=[.5, 1.08])
    plt.savefig(plotFile)
    plt.show()
    
# Function to reduce radial resolution for plot and tabulation
# Plot always uses 361 radials, tabular list shows maximum allowed
# for all radials surrounding the listed radial
# Input required - sorted list with required relative field a 1 degree
# resolution 0 through 360 degrees and desired radials - 72 or 36 
# Output: Sorted list with desired radials + 1 (0 repeated), dictionary with 361 values 
def ReducedPlotRadials(RelFieldList, radialCount):
    if radialCount == 72 :
        x = 0
        NewRadialDict = {}
        NewRadialList = []
        for NewAngle in range(0,359,5):
            RFtemp = []
            for x in range(0,6) :     # We overlap minimum search by one radial to possible problems with
                y = x - 3 + NewAngle  # IX on a fractional bearing due to TVStudy antenna pattern handlng
                if (y < 0) :
                    y = y + 360
                RFtemp.append(RelFieldList[y][1])
            NewRadialList.append(min(RFtemp))
#            print(NewAngle,min(RFtemp))
    elif radialCount == 36 :
        x = 0
        NewRadialDict = {}
        NewRadialList = []
        for NewAngle in range(0,359,10):
            RFtemp = []
            for x in range(0,12) :    # We overlap minimum search by one radial to possible problems with
                y = x - 5 + NewAngle  # IX on a fractional bearing due to TVStudy antenna pattern handlng
                if (y < 0) :
                    y = y + 360
                RFtemp.append(RelFieldList[y][1])
            NewRadialList.append(min(RFtemp))
#            print(NewAngle,min(RFtemp))
    else :
        NewRadialList = []
#        print("Radial reduction limited to 72 to 36 radials") 
    return(NewRadialList)
            
    
# Function to display 72 radial antenna azimuth patterns on a polar plot
def plotAzimuthPattern72(PatternList, plotFile) :
    angles = np.arange(0,361.0,5)
    PatternList.append(PatternList[0])
    Pattern1 = PatternList
#    for radial in angles :        # Put relative fields in a list start with 0 degree radial
#        Pattern1.append(PatternDict1[radial])        
    theta = (np.pi/180.0) * angles
    mpl.rc("ytick",labelsize = "small")
    fig1 = plt.figure(figsize = (5,5.4), dpi=120, facecolor='White')
    ax1 = fig1.add_axes([0.1,0.1,0.8,0.8], polar=True)
    ax1.set_theta_zero_location("N")
    ax1.set_theta_direction(-1)
    ax1.set_ylim(0,1)
    ax1.set_yticks(np.arange(0,1,0.1))
    ax1.plot(theta,Pattern1,lw=2.5)
    ax1.set_title(plotFile, position=[.5, 1.08])
    plt.savefig(plotFile)
    plt.show()

# Function to display 36 radial antenna azimuth patterns on a polar plot
def plotAzimuthPattern36(PatternList, plotFile) :
    angles = np.arange(0,361.0,10)
    PatternList.append(PatternList[0])
    Pattern1 = PatternList
#    for radial in angles :        # Put relative fields in a list start with 0 degree radial
#        Pattern1.append(PatternDict1[radial])        
    theta = (np.pi/180.0) * angles
    mpl.rc("ytick",labelsize = "small")
    fig1 = plt.figure(figsize = (5,5.4), dpi=120, facecolor='White')
    ax1 = fig1.add_axes([0.1,0.1,0.8,0.8], polar=True)
    ax1.set_theta_zero_location("N")
    ax1.set_theta_direction(-1)
    ax1.set_ylim(0,1)
    ax1.set_yticks(np.arange(0,1,0.1))
    ax1.plot(theta,Pattern1,lw=2.5)
    ax1.set_title(plotFile, position=[.5, 1.08])
    plt.savefig(plotFile)
    plt.show()

    
# Create list of interfering stations and power reduction required to protect them
# Create list with power reductions requires at defined intervals of 1, 5 and 10 degrees

results = getTVstudyData()
SourceInfoList = results[0]
sourceCall = SourceInfoList[2]
SourceFacID = SourceInfoList[6]
SourceFileNum = SourceInfoList[7]
SourceID = getSourceID(SourceFileNum)[0]
ListofVictims = results[1]
AllIxStudies = []

# ixCase [ixCall, ixStatus, ixScenarioNum, ExistingIxFreePop, StudyIxFreePop, ixService, ixFileNumber]
# filenumber = ixCase[6], scenarionumber = ixCase[2]
# Scenarios to study stored as dictionary with filenumber key and scenarios in a list
# Test: If (ixFilenum in StudyScenario) & (10 in StudyScenario[ixFilenum])


for i, ixCase in enumerate(ListofVictims):
    if (StudyScenario == {}) & (IgnoreScenario == {}):
        IxStudy = VictimIxStudy(SourceID, SourceInfoList,ixCase)
        AllIxStudies = AllIxStudies + IxStudy[4]
    elif (StudyScenario != {}) :
        ixFileNum = ixCase[6]
        ScenarioNum = ixCase[2]
        if ixFileNum in StudyScenario :
            if (StudyScenario[ixFileNum] == 'ALL') :
                IxStudy = VictimIxStudy(SourceID, SourceInfoList,ixCase)
                AllIxStudies = AllIxStudies + IxStudy[4]
            elif ((ixFileNum in StudyScenario) & (ScenarioNum in StudyScenario[ixFileNum])) :
                IxStudy = VictimIxStudy(SourceID, SourceInfoList,ixCase)
                AllIxStudies = AllIxStudies + IxStudy[4]
    elif (IgnoreScenario != {}) :
        ixFileNum = ixCase[6]
        ScenarioNum = ixCase[2]
        if ixFileNum in IgnoreScenario :
            if IgnoreScenario[ixFileNum] != 'ALL' :
                if not(ScenarioNum in IgnoreScenario[ixFileNum]) :
                    IxStudy = VictimIxStudy(SourceID, SourceInfoList,ixCase)
                    AllIxStudies = AllIxStudies + IxStudy[4]
        else :
            IxStudy = VictimIxStudy(SourceID, SourceInfoList,ixCase)
            AllIxStudies = AllIxStudies + IxStudy[4]
            
TxLatitude = SourceInfoList[3]
TxLongitude = SourceInfoList[4]
CombinedResults = CombineVictims(AllIxStudies,TxLatitude,TxLongitude)
CombinedRadials = CombinedResults[0]
maxReduction = CombinedResults[2]

CombinedCellsToClear = CombinedResults[1]
if (UseSimpleKML == True) :
    kml = simplekml.Kml(open=1)
    for Key, Point_Lat, Point_Lon, Point_Area, Point_Pop, Point_HH, U_Source_Key, U_Site_Num, D_Source_Key, D_Site_Num, D_Field_Str, Serv_Flag, U_Field_Str, Rcv_Pat_Adj, R2T_Bearing, DU_Ratio, DU_Thres, SNR_Fact, Causes_Ix, D_KWX_Flag, U_KWX_Flag, IX_level in CombinedCellsToClear :
        pnt = kml.newpoint()
        pnt.name = Point_Pop
        pnt.iconstyle.color = simplekml.Color.red
        pnt.description = 'Point_Pop: ' + Point_Pop + '   DU_Ratio: ' + DU_Ratio + '   DU_Thres: ' + DU_Thres + '   IX_level(dB): ' + '{:2.3f}'.format(float(IX_level)) + '   D_Source_Call: ' + sourcesDict[D_Source_Key][2] 
        pnt.coords = [('-'+Point_Lon, Point_Lat)]
    pnt = kml.newpoint(name=sourceCall)
    pnt.coords = [('-'+str(TxLongitude), str(TxLatitude))]
    pnt.style.labelstyle.color = simplekml.Color.yellow
    pnt.style.labelstyle.scale = 1.1
    kml.save("Cells_contributing_interference_over_deminimis_level_from "+sourceCall+".kml")

# Generate required antenna pattern file
RelFieldList = []
# print('SourceFacID for Pattern:  ',SourceFacID)
PatternFCC = getTVstudyPattern(SourceFacID)
PatternRequired = {}
Az360 = range(361)
if 0 in CombinedRadials :
    CombinedRadials[360] = CombinedRadials[0]
for radial in Az360 :
    if radial in CombinedRadials :
        PatternRequired[radial] = PatternFCC[radial]*math.pow(10,-(CombinedRadials[radial]/20))
    else :
        PatternRequired[radial] = PatternFCC[radial]
    RF = PatternRequired[radial]
    RelFieldList.append([radial,RF])
    RelFieldList = sorted(RelFieldList,key=itemgetter(0))

# Generate plot of required antenna pattern
if (UseMatplotlib == True):
    plotFile = sourceCall + "_Required_Pattern.png"
    plotAzimuthPattern(PatternRequired, plotFile)

# Create csv file suitable for import into TV Study
patternFile = sourceCall + "_Required_Pattern.csv"
ReqPatternCsv = open(patternFile,'w')
for item in RelFieldList :
    radial = '{:3.0f}'.format(float(item[0]))+','+'{:1.4f}'.format(float(item[1]))+'\n'
    ReqPatternCsv.write(radial)
ReqPatternCsv.close()


# Generate 72 radial antenna pattern
ReducedPlotList = ReducedPlotRadials(RelFieldList, 72)
# print(ReducedPlotList)

# Plot 72 radial antenna pattern
if (UseMatplotlib == True):
    plotFile = sourceCall + "_Required_Pattern-72-radials.png"
    plotAzimuthPattern72(ReducedPlotList, plotFile)

# Create 72 radial csv file suitable for import into TV Study
patternFile = sourceCall + "_Required_Pattern_72_radials.csv"
ReqPatternCsv = open(patternFile,'w')
for x in range(0,72) :
    radial = '{:3.0f}'.format(float(x * 5))+','+'{:1.4f}'.format(float(ReducedPlotList[x]))+'\n'
    ReqPatternCsv.write(radial)
ReqPatternCsv.close()

# Generate 36 radial antenna pattern
if (UseMatplotlib == True):
    ReducedPlotList = ReducedPlotRadials(RelFieldList, 36)
# print(ReducedPlotList)

# Plot 36 radial antenna pattern
if (UseMatplotlib == True):
    plotFile = sourceCall + "_Required_Pattern-36-radials.png"
    plotAzimuthPattern36(ReducedPlotList, plotFile)

# Create 36 radial csv file suitable for import into TV Study
patternFile = sourceCall + "_Required_Pattern_36_radials.csv"
ReqPatternCsv = open(patternFile,'w')
for x in range(0,36) :
    radial = '{:3.0f}'.format(float(x * 10))+','+'{:1.4f}'.format(float(ReducedPlotList[x]))+'\n'
    ReqPatternCsv.write(radial)
ReqPatternCsv.close()

ixstudy.write(ixstudyText)
ixstudy.close