# This is our primary file from which we analyze heat distribution in our system

import math
import numpy as np
#import fdsreader
#from fdsreader.bndf.utils import sort_patches_cartesian
import matplotlib.pyplot as plt
#import pandas as pd
import matplotlib
import matplotlib.animation as animation
import random

plt.style.use('_mpl-gallery-nogrid')

# for rectilinear model:
# Static position model defines grid as units of one with "adjacent"
# being defined as within a 2 by 2 cube centered on the particle in question 

# Hexagonal Spacial:
# Static positions defined as hexagonally offset
# This ensures the distance from each particle to 
# another is the same for all close particles

changeFactor = 1 #this is to make it easier to change the space between particles 
#This affects not only the placing but also the parameters for adjacency and 

class Particle(object): 
    # change __init__ structure to 3d 
    def __init__(self,xPos,yPos,zPos,temp): #these distances are in mm
        self.x = xPos
        self.y = yPos
        self.z = zPos
        self.t = temp
        
        self.Sx = 0
        self.Sy = 0
        self.Sz = 0
        
        #angular velocities
        self.wx = 0
        self.wy = 0
        self.wz = 0

        self.bondingGrade = -1
        
    def changeTemp(self,newTemp):
        self.t = newTemp
    
    def changeX(self,newX):
        self.x = newX
    
    def changeY(self,newY):
        self.y = newY
    
    def changeZ(self,newZ):
        self.z = newZ
        
    def changeVX(self,newXVel):
        self.Sx = newXVel
    
    def changeVY(self,newYVel):
        self.Sy = newYVel
        
    def changeVZ(self,newZVel):
        self.Sz = newZVel
    
    def changeWX(self,newXVel):
        self.wx = newXVel
    
    def changeWY(self,newYVel):
        self.wy = newYVel
        
    def changeWZ(self,newZVel):
        self.wz = newZVel
        
    def setGrade(self,grade):
        self.bondingGrade = grade
    
    def isPlate(self):
        return False


class Plate(Particle):
    def __init__(self,xPos,yPos,zPos): 
        self.x = xPos
        self.y = yPos
        self.z = zPos
        self.t = plateTemp
        
        #velocities
        self.Sx = 0
        self.Sy = 0
        self.Sz = 0
        
        #angular velocities
        self.wx = 0
        self.wy = 0
        self.wz = 0
        
        self.bondingGrade = 0

    #operates under the assumption of an unmoving plate of constant temperature
    def changeTemp(self,newTemp):
        self.t = plateTemp
        
    def changeVX(self,newXVel):
        self.Sx = 0
    
    def changeVY(self,newYVel):
        self.Sy = 0
        
    def changeVZ(self,newZVel):
        self.Sz = 0
        
    
    def changeWX(self,newXVel):
        self.wx = 0
    
    def changeWY(self,newYVel):
        self.wy = 0
        
    def changeWZ(self,newZVel):
        self.wz = 0
        
    def setGrade(self,grade):
        self.bondingGrade = 0
    
    def isPlate(self):
        return True


class Extruder(Particle): 
    
    pass 

class Material(Particle):
    
    pass 

#set up sim

#spacing is 1 unit for now

iLength = 4
kLength = 4
jLength = 2

npS = 3 #num per side 
basepS = 4 # num per side of base

#constriction (len = height)
heightChanges = np.array(["I","O","I"])
height = len(heightChanges) #doesn't include base

plateTemp = 50 #normal print bed temp
bondingTemp = 68.5 #glass point
extrusionTemp = 110 # C
rateOfEx = 3 #(particle per second)
timeFrame = 20 #seconds
frameNum = 100000
visc = 1287000 #viscosity N/m^2 (Pa)
Area = .000001 #area of side m^2 
thermalExpansion = 4.17*10**(-4) #% per C
TD = 5.8244*10**(-7)

speedCap = .02 # m/s so 20mm/s
forceCap = 100 #F  0.700348550209857

E = 17.5/0.672475 #base energy of bond (in Joules) such that the max force is 
m = .00000124 #kilograms
systemTime = 0.0 #for reference whenever an update is made to the system the time should be adjusted
repeating = False

'''
#just some particles
particleList = np.array([Particle(1,0,0,100),Particle(0,0,0,100),Particle(2,0,0,100),Particle(4,1,0,100),Particle(4,0,0,100)])
particleList[1].changeVY(.5)
'''
#Make a rectilinear set (rudimentary)
'''
particleList = np.array([None]*(iLength*kLength*jLength+(iLength+2)*(kLength+2)))
particleListHist = np.array([[None]*(iLength*kLength*jLength+(iLength+2)*(kLength+2))]*(frameNum+1))

i=0
while(i<iLength):
    k=0
    while(k<kLength):
        j=0
        while(j<jLength):
            realIndex = iLength*kLength*j+kLength*i+k
            #particleList[realIndex] = Particle(i+1,k+1,j+1,random.random()*100)
            particleList[realIndex] = Particle(i+1,k+1,j+1,0)
            #if(j==0):
                #particleList[realIndex].changeTemp(100)
            
            j+=1
        k+=1
    i+=1

i=0
while(i<iLength+2):
    k=0
    while(k<kLength+2):
        realIndex = (kLength+2)*i+k+iLength*kLength*jLength
        print(realIndex)
        particleList[realIndex] = Plate(i,k,0)
        k+=1
    i+=1
    

n=0
while(n<frameNum+1):
    i=0
    while(i<iLength):
        k=0
        while(k<kLength):
            j=0
            while(j<jLength):
                realIndex = iLength*kLength*j+kLength*i+k
                if(n==0):
                    particleListHist[n][realIndex] = Particle(i+1,k+1,j+1,particleList[realIndex].t)
                else:
                    particleListHist[n][realIndex] = Particle(i+1,k+1,j+1,None)
                j+=1
            k+=1
        i+=1
    i=0
    while(i<iLength+2):
        k=0
        while(k<kLength+2):
            realIndex = (kLength+2)*i+k+iLength*kLength*jLength
            particleListHist[n][realIndex] = Plate(i,k,0)
            k+=1
        i+=1    
    n+=1
'''

#Make a Hexagonal Set
p = True #if not in between
k = npS
setLen = 0
i=0
maxK = npS
#print("1 ", setLen)

while(i<len(heightChanges)):
    if(p):
        if(heightChanges[i]=="O"):
            k+=1
        #nothing happens to k on in from full
        setLen += 3*k**2-6*k+3
        p=False
        #print("1 ", setLen)
    else:
        if(heightChanges[i]=="I"):
            k-=1
        #nothing happens to k on out from half
        setLen += 3*k**2-3*k+1
        p=True
        #print("2 ", setLen)
    
    i+=1
    
setLen += 3*basepS**2-3*basepS+1
        
#print(setLen)

s2=math.sqrt(2)
s3=math.sqrt(3)

particleList = np.array([None]*setLen)
#particleListHist = np.array([[None]*setLen]*(frameNum+1))
ind = 0
l = 0
alternator = True
numIn = 0
heightFactorX = 0
heightFactorY = 0
while(l <= height):
    c = basepS
    npSs=basepS
    
    baseOffsetX = .5*(basepS-npS)
    baseOffsetY = s3/2*(basepS-npS)
    
    row = 0
    if(l==0): #for bottom plate
        while(row<2*npSs-1):
            i = 0
            while(i<c):
                if(row<npSs-1):
                    particleList[ind]=Plate(-1/2*row+i - baseOffsetX, s3/2*row-baseOffsetY, 0) #Need to orient properly
                    ind+=1
                else:
                    particleList[ind]=Plate(1-npSs+1/2*row+i - baseOffsetX, s3/2*row-baseOffsetY, 0)
                    ind+=1
                #print("Particle #: ",ind-1," x: ",particleList[ind-1].x, " y: ", particleList[ind-1].y, " z: ", particleList[ind-1].z)
                i+=1
            if(row<npSs-1):
                c+=1
            else:
                c-=1
            row+=1
    else:
        if(math.trunc(numIn/2)==numIn/2): #if numIn Even ie retains shape
            c=npS-numIn/2
            d=npS-1-numIn/2
            while(row<2*npS-1-numIn):
                i = 0
                while(i<c):
                    if(row<d):
                        particleList[ind]=Particle(-1/2*row+i + heightFactorX, s3/2*row+heightFactorY, l*s2/s3,extrusionTemp)
                        ind+=1
                    else:
                        particleList[ind]=Particle(1-npS+numIn/2+1/2*row+i + heightFactorX, s3/2*row+heightFactorY, l*s2/s3,extrusionTemp)
                        ind+=1
                    #print("Particle #: ",ind-1," x: ",particleList[ind-1].x, " y: ", particleList[ind-1].y, " z: ", particleList[ind-1].z, " row: ",row, " i: ", i, " c: ", c, " d: ", d)
                    i+=1
                    
                if(row<d):
                    c+=1
                else:
                    c-=1
                row+=1
        elif(not(alternator)): #should always be true since when entering into noneven numIn it should have forward structure
            c = npS-math.ceil(numIn/2)
            d = math.ceil(npS-(1+numIn)/2)
            while(row<2*npS-1-numIn):
                i = 0
                while(i<c):
                    if(row<d):
                        particleList[ind]=Particle(-1/2*row+i + heightFactorX, s3/2*row+heightFactorY, l*s2/s3,extrusionTemp)
                        ind+=1
                    else:
                        particleList[ind]=Particle(-1/2*2*d+1/2*row+i + heightFactorX, s3/2*row+heightFactorY, l*s2/s3,extrusionTemp)
                        ind+=1
                    #print("Particle #: ",ind-1," x: ",particleList[ind-1].x, " y: ", particleList[ind-1].y, " z: ", particleList[ind-1].z)
                    i+=1
                if(row<d):
                    c+=1
                else:
                    c-=1
                row+=1
        else:
            print("error: bad alternation")
    if(l<len(heightChanges)):
        if(heightChanges[l] == "I"):
            if(alternator):
                heightFactorX += 1/2
                heightFactorY += 1/s3/2
                alternator = False
                numIn +=1
            else:
                heightFactorY += 1/s3
                alternator = True
                numIn +=1
        else:
            if(alternator):
                heightFactorY -= 1/s3
                alternator = False
                numIn -=1
            else:
                heightFactorX -= 1/2
                heightFactorY -= 1/s3/2
                alternator = True
                numIn -=1
    l+=1
    


#make history file
'''
n=0
while(n<frameNum+1):
    i=0
    while(i<len(particleList)):
        if(n==0):
            if(i<3*n**2-3*n+1): #gives num of particles in base
                particleListHist[n][i]=Plate(particleList[i].x,particleList[i].y,particleList[i].z)
            else:
                particleListHist[n][i]=Particle(particleList[i].x,particleList[i].y,particleList[i].z,particleList[i].t)
        else:
            particleListHist[n][i]=Particle(None,None,None,None)
        i+=1
    n+=1
'''
#systemTimeList = np.array([0.0]))        

adjacentLists = np.array([[None]*30]*len(particleList)) #indexes
adjacentListsDists = np.array([[None]*30]*len(particleList)) #distances
adjacentListsPotential = np.array([[None]*30]*len(particleList)) #potentials
    
    


'''
i=0
while(i<len(particleList)):
    print("i: ", i)
    print("particle at x ",particleList[i].x, " y ", particleList[i].y, "with temp ", particleList[0].t)
    i+=1
'''


#Consider adjacent particles as defined by relevant area and thickness to transfer heat
#∆Q ~ Area of conduction / thickness of separating material (in this case distance between particles)* ∆T * t
#for particles that are getting farther away from one another, the area goes to zero and the 'thickness' goes to infinity
#think of the area at the percent of your view taken up by the particle and it getting smaller as it travels away 
    
#Since it is difficult to calculate and make sense of that kind of area in this particle based model
# There are two paths: Refining the grid while keeping well define conduction relations
# and considering a completely free interaction based on more accurate contact surfaces

#the first method
def runTimeStep1(systemTimeI):

    #Heat Exchange
    changeInHeat = np.array([0.0]*len(particleList))
    
    #Kinetic Exchange
    changeInVX = np.array([0.0]*len(particleList))
    changeInVY = np.array([0.0]*len(particleList))
    changeInVZ = np.array([0.0]*len(particleList))
    changeInWX = np.array([0.0]*len(particleList))
    changeInWY = np.array([0.0]*len(particleList))
    changeInWZ = np.array([0.0]*len(particleList))
    
    #adjacent is defined by with a distance less than the grid size (with additional margin)
    ind = -1
    k=0
    Fmax = 9.81*m 
    
    particleNum = math.floor(systemTimeI*rateOfEx)+3*basepS**2-3*basepS+1
    if(particleNum>len(particleList)):
        particleNum = len(particleList)
    while(k<particleNum):
        particle = particleList[k]
        lx = particle.x
        ly= particle.y
        lz = particle.z
        ind +=1
        #print("particle number ",k)
        #print("particles here ",particleNum)
        #for particle in particleList:
            #print("X: ",particle.x, " Y: ", particle.y, " Z: ", particle.z, " T: ", particle.t)
        
        #heat conduction into an infinite slab
        #if(lz == 1):
            #particle.changeTemp(100)

        #find adjacent
        
        #if bonded fully
        if(adjacentLists[ind][0]!=None):
            #print("bonded")
            i=0
            while(i<len(adjacentLists[0]) and adjacentLists[ind][i]!=None):
                particleActive = particleList[adjacentLists[ind][i]]
                dist = math.sqrt((particleActive.x-lx)**2+(particleActive.y-ly)**2+(particleActive.z-lz)**2)/1000 #in mm to convert to m
                heatOut = (particle.t-particleActive.t)/(dist**2)/2*TD
                changeInHeat[adjacentLists[ind][i]] += heatOut
                changeInHeat[ind] -= heatOut
                #increases general attraction
                #makes a new potential curve such that the bonded distance is the well and it's strength depends on where it is 
                dist = dist*1000 #back to mm
                F=12*E*((dist+1-adjacentListsDists[ind][i])**(-7)-(dist+1-adjacentListsDists[ind][i])**(-13))*4*adjacentListsPotential[ind][i] #4 is from baseline since E/4 gives the potential at the well
                if(F>Fmax):
                    Fmax=F
                changeInVX[adjacentIndices[i]] -= F/m*(particleActive.x-particle.x)/dist
                changeInVX[ind] += F/m*(particleActive.x-particle.x)/dist
                changeInVY[adjacentIndices[i]] -= F/m*(particleActive.y-particle.y)/dist
                changeInVY[ind] += F/m*(particleActive.y-particle.y)/dist
                changeInVZ[adjacentIndices[i]] -= F/m*(particleActive.z-particle.z)/dist
                changeInVZ[ind] += F/m*(particleActive.z-particle.z)/dist
                i+=1
                
                #Torque and Angular Momentum
                d = [(particleActive.x-particle.x)/1000,(particleActive.y-particle.y)/1000,(particleActive.z-particle.z)/1000] #distance vector
                v1 = [particleActive.Sx-particle.Sx,particleActive.Sy-particle.Sy,particleActive.Sz-particle.Sz] #bare velocity difference
                
                wa =[particleActive.wx,particleActive.wy,particleActive.wz] #angular velocities
                wh =[particle.wx,particle.wy,particle.wz]
                
                r1 = [d[0]/2,d[1]/2,d[2]/2] #'radius' vector of home particle
                r2 = [-r1[0],-r1[1],-r1[2]] #'radius' vector of active particle
                
                sva = np.cross(wa,r2) #surface velocity of active particle
                svh = np.cross(wh,r1) #surface velocity of home particle
                
                vtot =[v1[0]+sva[0]-svh[0],v1[1]+sva[1]-svh[1],v1[2]+sva[2]-svh[2]] #total velocity at the surface of the interaction from the reference of the active particle
                #print(vtot,"  Mag ",math.sqrt(vtot[0]**2+vtot[1]**2+vtot[2]**2))
                
                #Since friction coefficient is proportional to velocity difference
                friction=-visc/(math.sqrt(d[0]**2+d[1]**2+d[2]**2))*Area
                Fc = [vtot[0]*friction,vtot[1]*friction,vtot[2]*friction]
                FcM = math.sqrt(Fc[0]**2+Fc[1]**2+Fc[2]**2)
                if(FcM>Fmax):
                    Fmax=FcM
                
                
                Tq = np.cross(r1,Fc) #torque for home particle
                I = 2/5*m *(r1[0]**2+r1[1]**2+r1[2]**2) #I=2/5*mr^2
                #α = T/I 
                A = [Tq[0]/I,Tq[1]/I,Tq[2]/I]
                
                changeInWX[ind] += A[0]
                changeInWY[ind] += A[1]
                changeInWZ[ind] += A[2]
                
                #since the difference is in the orientation of the radius vector it's as simple as flipping a sign
                changeInWX[adjacentIndices[i]] -= A[0]
                changeInWY[adjacentIndices[i]] -= A[1]
                changeInWZ[adjacentIndices[i]] -= A[2]
                
                #Friction effect on the particles in general
                changeInVX[adjacentIndices[i]] -= Fc[0]
                changeInVX[ind] += Fc[0]
                changeInVY[adjacentIndices[i]] -= Fc[1]
                changeInVY[ind] += Fc[1]
                changeInVZ[adjacentIndices[i]] -= Fc[2]
                changeInVZ[ind] += Fc[2]
            if(particle.t > bondingTemp):
                #if not within temp requirements (too hot) and previously bonded wipe bonding
                #essentially the particle has remelted
                i=0
                while(i<len(adjacentLists[ind])):
                    adjacentLists[ind][i]=None
                    adjacentListsDists[ind][i]=None
                    adjacentListsPotential[ind][i]= None
                    particle.setGrade(-1)
                    i+=1
        else:   
            adjacentIndices = np.array([None]*30)
            adjDists = np.array([None]*30)
            nextEnter = 0
            i = 0
            while(i<particleNum):
                #general block to limit calculation complexity
                if((abs(particleList[i].x) - lx)<1.5 and (abs(particleList[i].y - ly)<1.5) and (abs(particleList[i].z-lz)<1.5) and not(ind == i)):
                    #print("here")
                    distanceTo = math.sqrt((particleList[i].x-lx)**2+(particleList[i].y-ly)**2+(particleList[i].z-lz)**2)
                    #print(distanceTo)
                    if(distanceTo < 1.3):
                        #print("here2")
                        adjacentIndices[nextEnter]=i
                        adjDists[nextEnter]=distanceTo
                        nextEnter +=1
                        if(nextEnter == 30):
                            print("Error: solid too dense -- too many 'adjacent' particles found")
                            i = len(particleList)

                i+=1
        
            i=0
            #print(nextEnter)
            while(i < nextEnter):
                particleActive = particleList[adjacentIndices[i]]
                particle = particleList[ind]
                dist = adjDists[i]/1000 # from mm to m
            
                #Temp
                #∆Q ~ 1cm^2 / math.sqrt((particleList[i].x-lx)**2+(particleList[i].y-ly)**2) * (tempLocal - tempExternal)* timeStep /2   Divided by two since all interactions will be considered twice from both sides
                heatOut = (particle.t-particleActive.t)/(dist**2)/2*TD # * some proportionality constant
            
                #change in still needs timeStep to be factored in in order to situationally lower it
                changeInHeat[adjacentIndices[i]] += heatOut
                changeInHeat[ind] -= heatOut
            
                #Kinetics
                #Force given from differentiated Lennard Jones Potential
                #F(d)=3E(x^(-7)-x^(-13))
           
                #We know impulse since we have a ∆t and Force
                #j = ∆p
                #m1v1+j=m2v2
                #v1+j/m=v2
                #j = F*∆t
                #∆v = F*∆t/m
                #Since force is vectorized, it can be rewritten as (with * * giving vector):
                # *v* = F*∆t/m * *r*/r
                #vx=F*∆t/m* ∆x/dist
                 
                TChange = particle.t/bondingTemp*4-3 #this for change with temperature
                dist = dist*1000 #back to mm
                F=12*E*(dist**(-7)-dist**(-13))/TChange
                #if(F>forceCap):
                    #F=forceCap
                #elif(F<-forceCap):
                    #F=-forceCap
                if(F>Fmax):
                    Fmax=F
                if(not(particleList[ind].isPlate())):
                    #print("index: ",ind, " active: ",adjacentIndices[i])
                    #print("dist ",dist)
                    #print("F ",F)
                    pass
                changeInVX[adjacentIndices[i]] -= F/m*(particleActive.x-particle.x)/dist
                changeInVX[ind] += F/m*(particleActive.x-particle.x)/dist
                changeInVY[adjacentIndices[i]] -= F/m*(particleActive.y-particle.y)/dist
                changeInVY[ind] += F/m*(particleActive.y-particle.y)/dist
                changeInVZ[adjacentIndices[i]] -= F/m*(particleActive.z-particle.z)/dist
                changeInVZ[ind] += F/m*(particleActive.z-particle.z)/dist
                
                
                #Torque and Angular Momentum
                d = [(particleActive.x-particle.x)/1000,(particleActive.y-particle.y)/1000,(particleActive.z-particle.z)/1000] #distance vector in mm converting to m
                v1 = [particleActive.Sx-particle.Sx,particleActive.Sy-particle.Sy,particleActive.Sz-particle.Sz] #bare velocity difference
                
                wa =[particleActive.wx,particleActive.wy,particleActive.wz] #angular velocities
                wh =[particle.wx,particle.wy,particle.wz]
                
                r1 = [d[0]/2,d[1]/2,d[2]/2] #'radius' vector of home particle
                r2 = [-r1[0],-r1[1],-r1[2]] #'radius' vector of active particle
                
                sva = np.cross(wa,r2) #surface velocity of active particle
                svh = np.cross(wh,r1) #surface velocity of home particle
                
                vtot =[v1[0]+sva[0]-svh[0],v1[1]+sva[1]-svh[1],v1[2]+sva[2]-svh[2]] #total velocity at the surface of the interaction from the reference of the active particle
                #print(vtot,"  Mag ",math.sqrt(vtot[0]**2+vtot[1]**2+vtot[2]**2))
                
                #drag equation + temp change for viscosity
                TChange = particle.t/bondingTemp*16-15
                friction=-visc/(math.sqrt(d[0]**2+d[1]**2+d[2]**2))/TChange*Area
                Fc = [vtot[0]*friction,vtot[1]*friction,vtot[2]*friction]
                FcM = math.sqrt(Fc[0]**2+Fc[1]**2+Fc[2]**2)
                if(FcM>Fmax):
                    Fmax=FcM
                
                Tq = np.cross(r1,Fc) #torque for home particle
                I = 2/5*m *(r1[0]**2+r1[1]**2+r1[2]**2) #I=2/5*mr^2
                #α = T/I 
                A= [Tq[0]/I,Tq[1]/I,Tq[2]/I]
                
                changeInWX[ind] += A[0]
                changeInWY[ind] += A[1]
                changeInWZ[ind] += A[2]
                
                #since the difference is in the orientation of the radius vector it's as simple as flipping a sign
                changeInWX[adjacentIndices[i]] -= A[0]
                changeInWY[adjacentIndices[i]] -= A[1]
                changeInWZ[adjacentIndices[i]] -= A[2]
                
                #Friction effect on the particles in general
                changeInVX[adjacentIndices[i]] -= Fc[0]
                changeInVX[ind] += Fc[0]
                changeInVY[adjacentIndices[i]] -= Fc[1]
                changeInVY[ind] += Fc[1]
                changeInVZ[adjacentIndices[i]] -= Fc[2]
                changeInVZ[ind] += Fc[2]
                
                i+=1
            
            #when a particle passes below a certain temperature it will hold on to its neighbors tightly but will not hold onto others tightly The tightness depends on the length of the bond
            if(particle.t < bondingTemp and not(particle.isPlate()) and nextEnter != 0): #glass point
                i=0
                while(i<nextEnter):
                    adjacentLists[ind][i]=adjacentIndices[i]
                    adjacentListsDists[ind][i]=adjDists[i]*(1-thermalExpansion*(extrusionTemp-bondingTemp))
                    adjacentListsPotential[ind][i]= (-adjDists[i]**(-6)/2+adjDists[i]**(-12)/4)
                    i+=1
                 
                particle.setGrade(gradeBonds(adjacentLists[ind],nextEnter,particle,adjacentListsPotential[ind]))
        k+=1
    i=0
    while(i<len(changeInVZ)):
        #∆v = F*∆t/m
        #F = mg
        #∆v = g*∆t
        #changeInVZ[i]-=9.81 #∆t held till later
        i+=1
    
    #set timeStep
    
    #find Vmax
    '''
    Vmax = 0
    i=0
    while(i< particleNum):
        speedSq = particleList[i].Sx**2+particleList[i].Sy**2+particleList[i].Sz**2
        if(speedSq>Vmax):
            Vmax = speedSq
        i+=1
    Vmax = math.sqrt(Vmax)
    '''
    timeStep = .000001
    #timeStep1 = .1*m/Fmax #find timestep such that the change in velocity is 100 mm/s
    #if(Vmax!=0):
        #timeStep2 = .0000001/Vmax #find timestep such that max distance moved is .001 mm --> enforces smooth consideration of potential
        #if(systemTime1>1):
            #print("V")
    #else:
        #timeStep2 = .0005
        #if(systemTime1>1):
            #print("F")
    #if(timeStep1<timeStep2):
        #timeStep = timeStep1
    #else:
        #timeStep = timeStep2
    #if(timeStep>.0005):
        #timeStep =.0005
    #if(not(particleList[i].isPlate())):
        #print("F ",Fmax)
        #print("Time ",systemTimeI)
    systemTime1 = timeStep + systemTimeI
    

    i = 0
    while(i< particleNum):
        #Heat
        particleList[i].changeTemp(particleList[i].t+changeInHeat[i]*timeStep)
        
        #Position
        if(not(particleList[i].isPlate())):
            #print("x pos ",particleList[i].Sx*timeStep+particleList[i].x)
            #print("y pos ",particleList[i].Sy*timeStep+particleList[i].y)
            #print("z pos ",particleList[i].Sz*timeStep+particleList[i].z)
            pass
        particleList[i].changeX(particleList[i].Sx*timeStep*1000+particleList[i].x) #from m/s to mm/s
        particleList[i].changeY(particleList[i].Sy*timeStep*1000+particleList[i].y)
        particleList[i].changeZ(particleList[i].Sz*timeStep*1000+particleList[i].z)
        
        
        #Velocity (deliberately calculated after position such that change in velocity is NOT incorporated into that timeStep's motion)
        speed = particleList[i].Sx+changeInVX[i]*timeStep
        #if(not(particleList[i].isPlate())):
            #print(speed)
        if(speed>speedCap): #Speed Caps to dampen oscillations
            particleList[i].changeVX(speedCap)
        elif(speed<-speedCap):
            particleList[i].changeVX(-speedCap)
        else:
            particleList[i].changeVX(particleList[i].Sx+changeInVX[i]*timeStep)
        speed = particleList[i].Sy+changeInVY[i]*timeStep
        #if(not(particleList[i].isPlate())):
            #print(speed)
        if(speed>speedCap): #Speed Caps to dampen oscillations
            particleList[i].changeVY(speedCap)
        elif(speed<-speedCap):
            particleList[i].changeVY(-speedCap)
        else:
            particleList[i].changeVY(particleList[i].Sy+changeInVY[i]*timeStep)
        speed = particleList[i].Sz+changeInVZ[i]*timeStep

        #if(not(particleList[i].isPlate())):
            #print(speed)
        if(speed>speedCap): #Speed Caps to dampen oscillations
            particleList[i].changeVZ(speedCap)
        elif(speed<-speedCap):
            particleList[i].changeVZ(-speedCap)
        else:
            particleList[i].changeVZ(particleList[i].Sz+changeInVZ[i]*timeStep)
        
        #Angular Momentum
        particleList[i].changeWX(particleList[i].wx+changeInWX[i]*timeStep)
        particleList[i].changeWY(particleList[i].wy*changeInWY[i]*timeStep)
        particleList[i].changeWZ(particleList[i].wz*changeInWZ[i]*timeStep)
        #The actual angular 'location' isn't important since the objects are rotationally symmetric
        
        
        i+=1
    
    return systemTime1


def gradeBonds(pList,nextEnter,home,pPotenList):
    score=0
    i=0
    while(i<nextEnter):
        if(not(particleList[pList[i]].isPlate())):
            score+=4*pPotenList[i]*50 # 50 points for good length bonds (4*potent has max 1)
            score+=(10-abs(particleList[pList[i]].t-bondingTemp))/2*5  # 50 for neighboring proximity to temp (yes can be negative but it shouldn't be)
        i+=1
    score = score/nextEnter #average
    return score


def plot():
    # plot
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

    x = np.array([None]*len(particleList))
    y = np.array([None]*len(particleList))
    z = np.array([None]*len(particleList))
    temps = np.array([None]*len(particleList))
    i=0
    while(i<len(particleList)):
        #print(i)
        x[i] = particleList[i].x
        y[i] = particleList[i].y
        z[i] = particleList[i].z
        temps[i] = particleList[i].t
        i+=1
    print("Scale Max: ",np.amax(temps)," Min: ", np.amin(temps))
    ax.scatter(x, y, z, c=temps, vmin=np.amin(temps), vmax = np.amax(temps))
    
    xStart = (-1-1/2*(npS-1))/changeFactor
    xEnd = (npS+1/2*(npS-1))/changeFactor
    yStart = -1/changeFactor
    yEnd = 2*npS/changeFactor
    
    print(xStart," - ", xEnd)
    print(yStart," - ", yEnd)

    #ax.set(xlim=(0, iLength+1), xticks=np.arange(0, iLength+1),
    #       ylim=(0, kLength+1), yticks=np.arange(0, kLength+1))
    #for haxagonal
    ax.set(xlim=(xStart, xEnd), xticks=np.arange(xStart, xEnd),
           ylim=(yStart, yEnd), yticks=np.arange(yStart, yEnd))

    plt.show()
    
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    
    print("Scale: -50 - 100")
    ax.scatter(x, y, z, c=temps, vmin=-50, vmax = 100)

    #ax.set(xlim=(0, iLength+1), xticks=np.arange(0, iLength+1),
    #       ylim=(0, kLength+1), yticks=np.arange(0, kLength+1),)
    #for hexagonal
    ax.set(xlim=(xStart, xEnd), xticks=np.arange(xStart, xEnd),
           ylim=(yStart, yEnd), yticks=np.arange(yStart, yEnd))

    plt.show()

#initial questioning and running algorithm
plot()
running = False
while(running):
    if(input("Would you like to run it for a timestep? y/n ")=="y"):
        runTimeStep1(systemTime)
        #systemTime += .1
        plot()
        '''elif(input("Would you like to run for a longer time? y/n ")=="y"):
        totalTime = (int)(input("How long would you like to run the simulation? (min)"))
        i = 0
        while(i < totalTime):
            runTimeStep1(.1)
            systemTime += .1
            i+=.1
        plot()
        '''
    else:
        running = False

#update particles to loop
def saveData(systemTime1):

        
    # remnants of animation
    #if repeating sequence (end of history is full)
    #if(not(particleListHist[frameNum][0].t==None)):
        #i=0
        #while(i<len(particleList)):
           # particleList[i] = particleListHist[num][i]
            #i+=1
    #else:
        #systemTimeList[num+1] = runTimeStep1(systemTimeList[num])
        #record
        #i=0
        
        #while(i<len(particleList)):
            #particleListHist[num+1][i].changeTemp(particleList[i].t)
            #particleListHist[num+1][i].changeX(particleList[i].x)
            #particleListHist[num+1][i].changeY(particleList[i].y)
            #particleListHist[num+1][i].changeZ(particleList[i].z)
            #i+=1
        #if(len(particleListHist)==num+2):
            #print out data
            #string1 = "currentWater%d.csv" % (int(percentGoal))
    particleNum = math.floor(systemTime*rateOfEx)+3*basepS**2-3*basepS+1
    if(particleNum>len(particleList)):
        particleNum = len(particleList)
    XHist=np.array([None]*particleNum)
    YHist=np.array([None]*particleNum)
    ZHist=np.array([None]*particleNum)
    THist=np.array([None]*particleNum)
    QHist=np.array([None]*particleNum)
    i=0
    while(i<particleNum):
        particleA = particleList[i]
        XHist[i]= particleA.x
        YHist[i]= particleA.y
        ZHist[i]= particleA.z
        THist[i]= particleA.t
        QHist[i]= particleA.bondingGrade  
        i+=1
    percentDone = math.trunc(systemTime/timeFrame*1000)/10
    string = 'dataX%f.csv'%(percentDone)
    np.savetxt(string,XHist, delimiter=",")
    string = 'dataY%f.csv'%(percentDone)
    np.savetxt(string,YHist, delimiter=",")
    string = 'dataZ%f.csv'%(percentDone)
    np.savetxt(string,ZHist, delimiter=",")
    string = 'dataT%f.csv'%(percentDone)
    np.savetxt(string,THist, delimiter=",")
    string = 'dataQ%f.csv'%(percentDone)
    np.savetxt(string,QHist, delimiter=",")
        
    #plot
    '''
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
    x = np.array([None]*particleNum)
    y = np.array([None]*particleNum)
    z = np.array([None]*particleNum)
    temps = np.array([None]*particleNum)
    i=0
    while(i<particleNum):
        #print(i)
        x[i] = particleList[i].x
        y[i] = particleList[i].y
        z[i] = particleList[i].z
        temps[i] = particleList[i].t
        i+=1
    ax.scatter(x, y, z, c=temps, vmin=0, vmax = 100)
    
    xStart = (-1-1/2*(npS-1))/changeFactor
    xEnd = (npS+1/2*(npS-1))/changeFactor
    yStart = -1/changeFactor
    yEnd = 2*npS/changeFactor
    zStart = 0 
    zEnd = height* .82
        
    ax.set(xlim=(xStart, xEnd), xticks=np.arange(xStart, xEnd),
            ylim=(yStart, yEnd), yticks=np.arange(yStart, yEnd),
            zlim=(zStart, zEnd), zticks=np.arange(zStart, zEnd))
    plt.show()
    '''
    return systemTime1
    '''
    # Some limit to see movement with static focus
    #zStart = -.1
    #zEnd = 1
    
    #ax.set(xlim=(0, iLength+1), xticks=np.arange(0, iLength+1),
    #       ylim=(0, kLength+1), yticks=np.arange(0, kLength+1))
    ''' # Z as said to be to see movement
  
        
    
    #systemTime = math.trunc(num*5+5)/frameNum # need to change to make time accurate

    #update plot  (remnants of plotting in animation during calculation)
    #print(systemTimeList[num])
    #print(systemTimeList[num+1])
    '''
    particleNum = math.floor(systemTime1*rateOfEx)+3*basepS**2-3*basepS+1
    if(particleNum>len(particleList)):
        particleNum = len(particleList)
    '''
    

numPrints = 5
start = 1 

#master loop
n=0
while(systemTime<timeFrame and n<frameNum):
    systemTime = runTimeStep1(systemTime)
    n+=1
    
    if(math.fmod(n,1000)==0):
        print("Completed ", n, " frames, time: ",systemTime)
        percentDone = math.trunc(systemTime/timeFrame*1000)/10
        print("Percent Done: ", percentDone)
    if(math.fmod(n,10000)==0):
        saveData(systemTime)
    if(timeFrame/systemTime>numPrints*start-1):
        saveData(systemTime)
        start+=1
        
#print(Fmax)
#print(Vmax)

length = len(particleList)
XHist=np.array([None]*length)
YHist=np.array([None]*length)
ZHist=np.array([None]*length)
THist=np.array([None]*length)
QHist=np.array([None]*length)
i=0
while(i<length):
    particleA = particleList[i]
    XHist[i]= particleA.x
    YHist[i]= particleA.y
    ZHist[i]= particleA.z
    THist[i]= particleA.t
    QHist[i]= particleA.bondingGrade  
    i+=1
    
string = 'dataXF%f.csv'%(systemTime)
np.savetxt(string,XHist, delimiter=",")
np.savetxt('dataYF.csv',YHist, delimiter=",")
np.savetxt('dataZF.csv',ZHist, delimiter=",")
np.savetxt('dataTF.csv',THist, delimiter=",")
np.savetxt('dataQF.csv',QHist, delimiter=",")
#np.savetxt('dataSystemTime.csv',systemTimeList, delimiter=",")

# plot once more with bond grades
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})

x = np.array([None]*len(particleList))
y = np.array([None]*len(particleList))
z = np.array([None]*len(particleList))
bg = np.array([None]*len(particleList))
i=0
while(i<len(particleList)):
    #print(i)
    x[i] = particleList[i].x
    y[i] = particleList[i].y
    z[i] = particleList[i].z
    bg[i] = particleList[i].bondingGrade
    i+=1
print("Scale Max: ",np.amax(bg)," Min: ", np.amin(bg))
ax.scatter(x, y, z, c=bg, vmin=np.amin(bg)-1, vmax = np.amax(bg)+2)
    
#xStart = (-1-1/2*(npS-1))/changeFactor
#xEnd = (npS+1/2*(npS-1))/changeFactor
#yStart = -1/changeFactor
#yEnd = 2*npS/changeFactor

xStart = -1
xEnd = 2
yStart = -3
yEnd = 4

    
print(xStart," - ", xEnd)
print(yStart," - ", yEnd)

#ax.set(xlim=(0, iLength+1), xticks=np.arange(0, iLength+1),
#       ylim=(0, kLength+1), yticks=np.arange(0, kLength+1))
#for haxagonal
ax.set(xlim=(xStart, xEnd), xticks=np.arange(xStart, xEnd),
        ylim=(yStart, yEnd), yticks=np.arange(yStart, yEnd))

plt.show()