##########################################################
### Import Necessary Modules

import argparse		           #provides options at the command line
import sys		               #take command line arguments and uses it in the script
import gzip		               #allows gzipped files to be read
import re		               #allows regular expressions to be used
import random                  #allows random numbers
import math                    #allows common math functions to be used

##########################################################
### Command-line Arguments
parser = argparse.ArgumentParser(description="A script to model virtual particles with particles as they pass through a double slit.")
parser.add_argument("-size", help = "The particle size, default = 1", default=1)
parser.add_argument("-slit", help = "The slit size (multiple of particle size, inf. height), default = 10", default=10)
parser.add_argument("-num", help = "The number of particles, default = 1000", default=1000)
parser.add_argument("-ang", help = "The maximum angle of virtual particle pairs (only allows 1-90), default = 90", default=90)
parser.add_argument("-mang", help = "The minimum angle of virtual particle pairs (only allows 1-90), default = 1", default=0)
parser.add_argument("-inter", help = "The size of the inter slit distance (multiple of particle size, inf. height), default = 5", default=5)
args = parser.parse_args()


class emitParticle():

    def __init__ (self):
        """Determine the initial trajectory of each particle (only one slit).  Assumes they are from one plane."""
        self.position = random.uniform(0.5*(float(args.size)),float(args.slit) - 0.5*(float(args.size)))
        self.angle = random.uniform(0,90)
        self.virtualParticle(self.position, self.angle)

    def virtualParticle(self, pos, ang):
        """Interaction of particle during transit through slit and after with the double slit"""
        #print("{}\t{}".format(pos, ang)) ### initial position and angle
        ### Left most slit wall position = 0
        if float(ang) < 90:  ###Assuming that the other angles will have similar patterns
            ### since this is a right triangle, we can solve using hypotenuse = adjacent/cos(angle)
            self.distanceL = float(pos)/math.cos(math.radians(ang)) ### The position is halfway into the particle
            self.distanceR = (float(args.slit) - float(pos))/math.cos(math.radians(ang))

            ### distance in virtual particles is just the number of pairs that can fit on the lines
            self.distanceLvpf, self.distanceLvp = math.modf((float(self.distanceL) - 0.5*float(args.size))/(float(args.size)*2)) ### need to include the size of the particle (left = 0.5)
            self.distanceRvpf, self.distanceRvp = math.modf((float(self.distanceR) - 0.5*float(args.size))/(float(args.size)*2)) ### need to include the size of the particle (right = 0.5)
            #print("{}\t{}\t{}\t{}".format(self.distanceLvp, self.distanceLvpf, self.distanceRvp, self.distanceRvpf))

            ### Calculate the posible change in distance (assuming that as long as the slit distance is small, the change should not be dependent on the distance but the virtual particle number fit)  
            ### need to find the new position for left and right (adjacent) and then mid of the two new distances
            ### adjacent = hypotenuse*cos(angle)
            ###Left side
            self.newHypotenuseL = (2*float(args.size)*(float(self.distanceLvp) + 1)) + 0.5*float(args.size)
            if float(self.distanceLvpf) == 0:
                self.newHypotenuseL = (2*float(args.size)*(float(self.distanceLvp))) + 0.5*float(args.size)
            self.newAdjacentL = float(self.newHypotenuseL)*math.cos(math.radians(ang))

            ###Right side
            self.newHypotenuseR = (2*float(args.size)*(float(self.distanceRvp) + 1)) + 0.5*float(args.size)
            if float(self.distanceRvpf) == 0:
                self.newHypotenuseR = (2*float(args.size)*(float(self.distanceRvp))) + 0.5*float(args.size)
            self.newAdjacentR = float(self.newHypotenuseR)*math.cos(math.radians(ang))
            #print("\t{}\t{}".format(self.newAdjacentL, self.newAdjacentR))

            ###Calculate the change in position from each side
            self.deltaL = (float(self.newAdjacentL) - float(pos)) ###Should be a positive number
            self.deltaR = float(self.newAdjacentR) - (float(args.slit)-float(pos)) ###Should be a positive number

            ###Calculate the distance in the center of the changes in positions
            self.finalPos1 = float(pos) + (float(self.deltaL) - float(self.deltaR))/2  ###Not exactly sure how to model this, part of the change in length is incorporated into the size of the deltas, so 
            ###half of the difference might make sense


            #####################################################################################
            #### Second slit influence ##########################################################
            #####################################################################################

            ###Left side
            self.adjSlitL = float(args.slit) - float(pos) + float(args.inter)
            self.distanceL2 = float(self.adjSlitL)/math.cos(math.radians(ang))
            self.distanceL2vpf, self.distanceL2vp = math.modf(float(self.distanceL2)/(float(args.size)*2))
            ### This needs to adjust for the other slit
            self.newHypotenuseL2 = (2*float(args.size)*(float(self.distanceL2vp) + 1)) + 0.5*float(args.size)
            if float(self.distanceL2vpf) == 0:
                self.newHypotenuseL2 = (2*float(args.size)*(float(self.distanceL2vp))) + 0.5*float(args.size)
            self.newAdjacentL2 = float(self.newHypotenuseL2)*math.cos(math.radians(ang))
            self.deltaL2 = float(self.newAdjacentL2) - float(self.adjSlitL) ###Should be a positive number

            ### Right side       
            self.adjSlitR = float(args.slit) - float(pos) + float(args.inter) + float(args.slit)
            self.distanceR2 = float(self.adjSlitR)/math.cos(math.radians(ang))
            self.distanceR2vpf, self.distanceR2vp = math.modf(float(self.distanceR2)/(float(args.size)*2))
            ### This needs to adjust for the other slit
            self.newHypotenuseR2 = (2*float(args.size)*(float(self.distanceR2vp) + 1)) + 0.5*float(args.size)
            if float(self.distanceR2vpf) == 0:
                self.newHypotenuseR2 = (2*float(args.size)*(float(self.distanceR2vp))) + 0.5*float(args.size)
            self.newAdjacentR2 = float(self.newHypotenuseR2)*math.cos(math.radians(ang))
            self.deltaR2 = float(self.newAdjacentR2) - float(self.adjSlitR) ###Should be a positive number

            ###Calculate the final pos
            self.finalPos2 = float(pos) + (float(self.deltaL) - float(self.deltaR))/2 + (float(self.deltaL2) - float(self.deltaR2))/2 ###Not sure how to model this properly

            #print("{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}".format(pos, ang, self.distanceL, self.distanceR, self.distanceLvp, self.distanceRvp, self.distanceLvpf, self.distanceRvpf, self.newAdjacentL, self.newAdjacentR))
            print("{}\t{}\t{}\t{}".format(pos, ang, self.finalPos1, self.finalPos2))
            

if __name__ == '__main__':
    print("{}\t{}\t{}\t{}".format("Initial.Positions","angle","finalPos1","finalPos2"))
    for i in range(0, int(args.num), 1):
        open_map = emitParticle()