######################### Stirling engine Code  #######################################
#Nikolai Lesack Feb/15/2022


# This code can be used to calculate the area inside a set of points
#To run this code you will want to modify the inputted data file, the data columns used 
#and the number of rectangle bins used termed "bins"

import numpy as np
import matplotlib.pyplot as plt

#set figure size
plt.rcParams['figure.figsize'] = [10, 10]

################ Import data ###############################################################
#Here are a few example data files that can be run by uncommenting them 
#To run your own code modify the file called

#Data = np.loadtxt("triangle.txt")
#Data = np.loadtxt("ellipse_noisy.txt")
Data = np.loadtxt("circle.txt")
#Data = np.loadtxt("noFrictionPV.dat")

#Will have to check which column is x and y when importing 
#data and change the x and y indices appropriately 
np.shape(Data)
x = Data[:,0]
y = Data[:,1]


#plot data 
plt.scatter(x, y)
plt.xlabel('$\Delta V$ (cm$^3$)', fontsize=20)
plt.ylabel('$\Delta P$ (kPa)', fontsize=20)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)


################### Calculating the area ##############################################

#this specifies the number of rectangels that will be used in the calculation, 
#the greater the number of rectangles the higher the accuracy however
#if the rectangles become too small an error will occur due to no points being 
#inside the rectangle bins
bins = 35 

#Next want to find the length of the area in the x-axis
#Find the largest x value 
min_x = min(x) 
max_x = max(x) 

#Next define the step size for each bin based length of shape in x and number of bins
x_step = (max_x - min_x)/bins 

#preset area for iterative calc
area = 0 



#Next use for loop to calculate area of each box
for l in range(bins):
    
    #set boundaries of the given rectangle bin
    x1 = min_x + x_step*l
    x2 = min_x + (l+1)*x_step
    
    #define empty list that will hold y-values in the bin
    y_list = [] 
    
    #This loop checks if a given value of x is in the current rectangular bin
    #and if it is it adds the corresponding y value to a list
    for q in range(len(x)): 
        
        if x1 <= x[q] < x2:
            y_list.append(y[q]) 
        
        
    #Next for the points in the given bin, find the mean y value
    y_mean = sum(y_list)/len(y_list)  
    

    #Seperate y values into those above and below range
    yhigh = [i for i in y_list if i >= y_mean] 
    ylow =  [i for i in y_list if i < y_mean]
 

    #Find mean of sperated y values 
    yhigh_mean = sum(yhigh)/len(yhigh)
    ylow_mean = sum(ylow)/len(ylow)
    
    #calculate area of rectangle and add it to running average 
    area = area + ((yhigh_mean - ylow_mean) * x_step)
    
    
    #Plot the rectanlge bin
    plt.plot([x1, x1], [yhigh_mean, ylow_mean], 'k--', lw=1)
    plt.plot([x2, x2], [yhigh_mean, ylow_mean], 'k--', lw=1)
    plt.plot([x1, x2], [yhigh_mean, yhigh_mean], 'k--', lw=1)
    plt.plot([x1, x2], [ylow_mean, ylow_mean], 'k--', lw=1)
    
plt.savefig('Jupyter PV area.pdf', format='pdf')
    
print('area =', area)

plt.show()