# ECG 
#
from visual import *
from visual.graph import *  # import graphing features
#

############################################
#
#       Student entered parameters

RI  = 0.62       # Enter the distance RI in meters as defined in the lab writeup here
RII = 1.09       # Enter the distance RII in meters as defined in the lab writeup here
FILENAME = 'SAK10sLeadIandII.csv'
#
#############################################
# Set this parameter to one if you want a curve of the
# dipole moment trajectory, set to 0 otherwise
#
plotcurve = 0

#############################################
#
# timestep in seconds
dt = 0.001
#
#############################################
#
# read the data
#
x=[]
f = open(FILENAME, 'r')
for ii in range(10000):
    x.append(f.readline())
#
#############################################
#
# parse the data
#
#create an array to hold the data
#
data = []
#
for ii in range(2,len(x)):
    #each element of x is a string, three comma separated values
    #ci is the comma index, a list of two positions in the
    #string where the commas are
    ci=[]
    for jj in range(len(x[ii])-1):
        s = x[ii][jj]
        if s==',':
            ci.append(jj)
    s2 = x[ii]

    if len(ci)>0:
        triple =[]
        triple.append(float(s2[0:ci[0]]))
        triple.append(float(s2[ci[0]+1:ci[1]]))
        triple.append(float(s2[ci[1]+1:len(s2)-2]))

        data.append(triple)
    
print "read data"
print shape(data)



#Now plot the data - we need two objects, a dipole arrow and
#a plot of the voltages from the two leads

# For plotting, the wireframe dimensions along Cartesian axes
oxa=20                         
oya=20
oza=20
# Initialize graphics scene - put the 3D box at x=50, y=0, make it 800 wide and 200 high
scene=display(width=500,height=500,x=0, y=0)

scene.ambient=0.24
scene.background=color.white
scene.title="Dipole moment of heart"
scene.center=(0,0,0)
scene.range=(oxa,oya,oza)
#this sets the viewing angle
scene.forward=(0,0,-1)


# set up the gdisplay options, with width = 1200 and height=600 and x=0 and y=500
#the graph basically covers the bottom half of the screen
graph1 = gdisplay(x=50, y=700, width=1200, height=300,
          title='Voltage vs. time: LEAD I ', xtitle='t (s)', ytitle='Potential (Volts)',
          foreground=color.black, background=color.white)

#this is a trick - vpython autoscales so plot white (invisible) data first, then plot
#actual data
functa = gcurve(gdisplay=graph1,delta=0.05, color=color.white)
funct1 = gcurve(gdisplay=graph1,color=color.black) 


graph2 = gdisplay(x=50, y=400, width=1200, height=300,
          title='Voltage vs. time: LEAD II ', xtitle='t (s)', ytitle='Potential (Volts)',
          foreground=color.black, background=color.white)

functb = gcurve(gdisplay=graph2,color=color.white) 
funct2 = gcurve(gdisplay=graph2,delta=0.05, color=color.red)

#this is a trick - vpython autoscales so plot white (invisible) data first, then plot
#actual data
for n in range(len(data)):

    functa.plot(pos=(data[n][0],data[n][1]))   
    functb.plot(pos=(data[n][0],data[n][2]))

#now define an arrow
dipole = arrow(pos=(0,0,0), axis = (1,1,0),shaftwidth=1,color=color.red)
refx = arrow(pos=(0,0,0), axis = (5,0,0),shaftwidth=0.1,color=color.blue)
refy = arrow(pos=(0,0,0), axis = (0,-5,0),shaftwidth=0.1,color=color.blue)
#define how many points to smooth over
smooth = 20
#make a list to compute the running average
runavgx = []
runavgy = []

####
# Make a heart shaped curve
#
# x= sin(t)cos(t)ln(|t|)
# y = sqrt(|t|)cos(t)
# -1<t<1
#
hearty = []
points=200
heartscale=20.0
origin = [0.0,-5.0]
for ii in range(points):
    t = -1.0 + 2.1*float(ii)/float(points)
    #print t,heartscale,sin(t),cos(t),log(abs(t)),origin[0],scale,sqrt(abs(t)),cos(t),origin[1]
    coords = [heartscale*sin(t)*cos(t)*log(abs(t))+origin[0],heartscale*sqrt(abs(t))*cos(t)+origin[1],0]
    hearty.append(coords)

curve(pos=hearty,color=color.red)

trace=[]
for n in range(len(data)):
    #make the loop execute 10 times per second - a 100 times slowdown
    rate(200)
    #compute the dipole moment
    #px = 0.001*data[n][0]*80.0*RI*RI
    #py = 0.001*data[n][1]*80.0*RII*RII/(sqrt(3.0)) - px/sqrt(3.0)    

    px = data[n][1]*80.0*RI*RI
    py = -(data[n][2]*80.0*RII*RII/(sqrt(3.0)) - px/sqrt(3.0))    

    
    runavgx.append(px)
    runavgy.append(py)

    if len(runavgx)>smooth:
        runavgx.remove(runavgx[0])
        runavgy.remove(runavgy[0])

    pxavg = sum(runavgx)/float(len(runavgx))
    pyavg = sum(runavgy)/float(len(runavgy))

    dipole.axis = (pxavg,pyavg,0)

    if plotcurve == 1:
        trace.append([pxavg,pyavg])
        curve(pos=trace,color=color.green)

    #plot the result of lead I measurement
    funct1.plot(pos=(data[n][0],data[n][1]))
    #plot the result of lead II measurement
    funct2.plot(pos=(data[n][0],data[n][2]))