python_postprocessing.py

import numpy as N
import re

def proces_raw_results(rawfile, savedir):

    '''
    rawfile - the 'simul' file that generated by SOLSTICE
    savedir - the directory for saving the organised results
    '''
    with open(rawfile) as f:
        content = f.read().splitlines()
        # content is a list with lots of string
    f.close()

    results=N.array([])

    # sun direction
    sun=re.findall("[-+]?\d*\.\d+|\d+", content[0])
    azimuth=sun[0]
    elevation=sun[1]

    # Global results
    line=list(map(float, content[1].split(' ')))
    num_hst=line[2]
    num_rays=line[3]

    potential=list(map(float, content[2].split(' ')))[0]#W
    potential_err=list(map(float, content[2].split(' ')))[1]

    absorbed=list(map(float, content[3].split(' ')))[0]
    absorbed_err=list(map(float, content[3].split(' ')))[1]    

    Fcos=list(map(float, content[4].split(' ')))[0]
    Fcos_err=list(map(float, content[4].split(' ')))[1] 

    shadow_loss=list(map(float, content[5].split(' ')))[0]
    shadow_err=list(map(float, content[5].split(' ')))[1]  
  
    missing_loss=list(map(float, content[6].split(' ')))[0]
    missing_err=list(map(float, content[6].split(' ')))[1] 

    material_loss=list(map(float, content[7].split(' ')))[0]
    material_err=list(map(float, content[7].split(' ')))[1]   
   
    atmospheric_loss=list(map(float, content[8].split(' ')))[0]
    atmospheric_err=list(map(float, content[8].split(' ')))[1]  


    # Target (receiver)
    target=content[9].split()
    rec_area=float(target[2]) # m2  
    rec_income=float(target[3])
    rec_income_err=float(target[4])
    
    rec_absorbed=float(target[13])
    rec_absorbed_err=float(target[14])
    rec_eff=float(target[23])
    rec_eff_err=float(target[24])


    #Virtual target
    virtual=content[10].split()
    vir_area=float(virtual[2])
    vir_income=float(virtual[25]) # back face available
    #print(vir_income)
    vir_income_err=float(virtual[4])
    
    Qtotal=potential
    Qcos=Qtotal*(1.-Fcos)
    Qshade=shadow_loss
    Qfield_abs=(Qtotal-Qcos-Qshade)*(1.-0.9)  # the mirror absorptivity needs to be adjusted here!!!!
    Q_att=atmospheric_loss
    Qabs=absorbed
    Qspil=vir_income-Qabs
    Qrefl=rec_income-Qabs
    Qblock=Qtotal-Qcos-Qshade-Qfield_abs-Qspil-Qabs-Qrefl-Q_att
    
    Q=N.array([Qtotal,Qcos,Qshade,Qfield_abs,Qblock,Q_att,Qspil,Qrefl,Qabs])
    #print(Q)
    efficiency_total=Qabs/Qtotal
    efficiency_exclu_inter=(Qabs+Qspil)/Qtotal
   # print('Interception: '+str(Qabs/(Qspil+Qabs))
    return efficiency_total,Q,efficiency_exclu_inter

def get_heliostat_to_receiver_data(simul, DNI, receiver_name):
	'''
	Reads all the information necessary to characterise heliostat_wise results from the simul file.
	Output is an array in which each line is the same heliostat as in the ioriginal .csv field file and the columns are the following information:
	heliostat_index, cosine_efficiency, shading loss (W), incident radiation on the heliostat(w), power delivered to the receiver (but not absorbed by the receiver) (w), heliostat absorption loss (W), atmospheric attenuation loss (W)
	'''

	# Script to exctract the efficiency and power from the helisotats to the receiver.
	simul_data = open(simul,'r')
	data = simul_data.readlines()
	simul_data.close()
	counts = [int(x) for x in data[1].split(' ')]
	
	# Find the receiver index:
	ridxs = 2+counts[0], 2+counts[0]+counts[1]

	recs = N.array([' '.join(x.split(' ')).split() for x in data[ridxs[0]:ridxs[1]]])
	
	recID = int(recs[recs[:,0]==receiver_name,1][0])
	

	# Get cosine efficiency and shadowing loss from Primaries data:
	pidxs = N.sum(counts[:2])+2, N.sum(counts[:3])+2
	hstats = N.array([' '.join(x.split(' ')[1:]).split() for x in data[pidxs[0]:pidxs[1]]], dtype=float)

	# Find the original line on the original field .csv file by processing the name:
	h_idx_orig = [int((x.split('.')[0])[2:]) for x in data[pidxs[0]:pidxs[1]]]

	
	areas = hstats[:,1]
	cosine_eff = hstats[:,3]
	q_shad = hstats[:,5]
	q_sun = DNI*areas
	q_in_h = DNI*areas*cosine_eff-q_shad

	# Get energy to receiver data from receiverXprimary data:
	rXpidxs = pidxs[1], pidxs[1]+counts[1]*counts[2]
	rXp = N.array([' '.join(x.split(' ')).split() for x in data[rXpidxs[0]:rXpidxs[1]]], dtype=float)
	r_mask = rXp[:,0]==recID
	#print(r_mask)
	
	# front face
	q_in_r = rXp[r_mask,2]
	q_abs_h = rXp[r_mask,8]
	q_atm_h = rXp[r_mask,10]
	'''
	# back face
	q_in_r = rXp[r_mask,22]
	q_abs_h = rXp[r_mask,28]
	q_atm_h = rXp[r_mask,30]
'''
	results = N.zeros((counts[2], 7))
	results[:,0] = h_idx_orig
	
	results[:,1] = cosine_eff
	results[:,2] = q_sun
	results[:,3] = q_in_h
	results[:,4] = q_in_r
	results[:,5] = q_abs_h
	#print(results[:,0])
	results = results[results[:,0].argsort()]
	return results[:,2],results[:,4]
	

def determine_field(simulfile, DNI, receiver_name, target_aperture_power, csv, csv_new, option=None):
	'''
	Function to select heliostats from a field based on an efficiency metric. The helisotat to receiver data are read and then enough heliostats are selected to match the input power required at the receiver. For now there are two options:
	- Cosine: selects heliostats based on their cosine efficiency
	- None: standard option to select heliostats based on their overall efficiency
	Returns: 
	The selected heliostat field data in the same format as the get_heliostat_to_receiver_data() function and ranked by efficiency.
	'''
	helio_data = get_heliostat_to_receiver_data(simulfile, DNI, receiver_name)
	
	if option == 'cosine':
		helio_eff = helio_data[:,1]
	else:
		helio_eff = helio_data[:,4]/helio_data[:,2]
	helio_data_ranked = helio_data[N.argsort(helio_eff)[::-1]]
	helio_eff=-N.sort(-helio_eff) # sort the hst eff

	cumul_p = N.add.accumulate(helio_data_ranked[:,4])
	helio_data_partial_ranked = helio_data_ranked[~(cumul_p>target_aperture_power)]
	# write all the chosen index 
	new_idx_org=helio_data_partial_ranked[:,0]
	num_hst_new=len(new_idx_org)
	
	# Script to exctract the information from the old csv file
	hst_info=N.loadtxt(csv,delimiter=',', skiprows=2)
	num_hst=len(hst_info)
	idx_array=N.arange(0,num_hst)
	hst_info_index=N.zeros((num_hst,8)) # add hst index to the first column of the array
	hst_info_index[:,0]=idx_array
	hst_info_index[:,1:]=hst_info
	left_hst = N.in1d(hst_info_index[:,0], new_idx_org)
	
	hst_info_new=hst_info_index[left_hst,:]

	pos_and_aiming=hst_info_new[:,1:]
	title=N.array(['x', 'y', 'z', 'foc', 'aim x', 'aim y', 'aim z', 'm', 'm', 'm', 'm', 'm', 'm', 'm'])
	pos_and_aiming=N.append(title, pos_and_aiming)
	pos_and_aiming=pos_and_aiming.reshape(len(pos_and_aiming)/7, 7)
	# plot
	plot(DM=15.70,pos_and_aiming=hst_info)
	plot(DM=15.70,pos_and_aiming=pos_and_aiming)
	N.savetxt(csv_new, pos_and_aiming, fmt='%s', delimiter=',')

def plot(DM, pos_and_aiming):
	num_hst_new=len(pos_and_aiming)-2
	print(num_hst_new)
	# plot 
	fig, ax = plt.subplots(figsize=(16,9))
	
	for i in range(2,num_hst_new):
		#print(pos_and_aiming[i,0], pos_and_aiming[i,1])
		circle = plt.Circle((float(pos_and_aiming[i,0]), float(pos_and_aiming[i,1])), radius=DM/2., ec='b',linewidth=0.5)
		plt.gca().add_patch(circle)
	plt.grid(linestyle='--')
	plt.xlim(-500.,500.)
	plt.ylim(00.,1000.)
	ax.tick_params(axis='both', which='major', labelsize=14)
	ax.set_aspect(1)
	plt.show()
	

	

if __name__=='__main__':
    simul='/home/admin-shuang/Desktop/SOLSTICE_relevant/Field_design/test_case/vtk/simul'
    csv='/home/admin-shuang/Desktop/SOLSTICE_relevant/Field_design/test_case/pos_and_aiming.csv'
    csv_new = '/home/admin-shuang/Desktop/SOLSTICE_relevant/Field_design/test_case/pos_and_aiming2.csv'
    DNI=900. #W/m2
    receiver_name='plate'
    target_aperture_power=60.e6 #W
    determine_field(simul, DNI, receiver_name,target_aperture_power,csv,csv_new)