import numpy as np 
import matplotlib.pyplot as plt 
import scipy
plt.rcParams.update({'font.size': 22})

data=np.loadtxt('PR2SpotVsCircPower5July2025.txt')
Expt_PR2_Spots=data[0,:]
Expt_CIRC_Powers=data[1,:]


data_19Feb=np.loadtxt('CleanExptData_19_Feb_2025.txt')


Pr2_Spot_19Feb=data_19Feb[0,:]
Power_19Feb=data_19Feb[1,:]
# print(data_19Feb.shape)



def setGeom(dtheta1):
	global X_prm,Y_prm, X_pr3, Y_pr3 , l1, l2, theta1_o, theta2_o, X_topEdge,Y_topEdge, X_bottomEdge, Y_bottomEdge
	l1= 15000
	l2= l1




	theta1_o=(0.339)*(np.pi/180)
	theta2_o=(1.264)*(np.pi/180)

	# dtheta1=-(0)*(np.pi/180)
	# dtheta=(.1)*(np.pi/180)

	X_prm=-l1*np.cos(theta1_o)
	Y_prm=-l1*np.sin(theta1_o)


	X_pr3=-l2*np.cos(theta2_o)
	Y_pr3=l2*np.sin(theta2_o)



	# BAFFLE Edges (TOP EDGE) 
	X_topEdge=-484.88 + (2.125*25.4)
	Y_topEdge=1.18 + (1.38*25.4)

	X_bottomEdge=-484.88 - (2.125*25.4)
	Y_bottomEdge=1.18 - (1.38*25.4)
	plt.figure(1)
	plt.plot([X_topEdge,X_bottomEdge],[Y_topEdge,Y_bottomEdge],marker='o',color='c')
	plt.plot([0,0],[-30,30],color='k',linewidth=2)
	offset=15
	plt.text(X_prm+offset, Y_prm+offset, 'PRM',size=15)
	plt.text(X_pr3-offset, Y_pr3-offset, 'PR3',size=15)
	plt.text(0+3*offset, 0-offset, 'PR2',size=15)
	# plt.arrow(X_topEdge,Y_topEdge,10, 10,width= 4) 
	plt.xlabel('X axis (mm)')
	plt.ylabel('Y axis (mm)')



def findIntersection(X1,Y1,X2,Y2,X3,Y3,X4,Y4):
	xTop_Det=np.array([[X1*Y2-X2*Y1,X1-X2],[X3*Y4-X4*Y3,X3-X4]])
	Bot_Det=np.array([[Y2-Y1,X1-X2],[Y4-Y3,X3-X4]])

	yTop_Det=np.array([[Y2-Y1,X1*Y2-X2*Y1],[Y4-Y3,X3*Y4-X4*Y3]])

	X_int=np.linalg.det(xTop_Det)/np.linalg.det(Bot_Det)
	Y_int=np.linalg.det(yTop_Det)/np.linalg.det(Bot_Det)
	return X_int,Y_int





def findNormalDistance(X0,Y0,X1,Y1,X2,Y2):

	Den=np.sqrt((X1-X2)**2+(Y2-Y1)**2)
	Num=(Y2-Y1)*X0+(X1-X2)*Y0+X2*Y1-Y2*X1
	X_Foot= ( (Y2-Y1)*(X1*(Y2-Y0)+X2*(Y0-Y1)) + X0*((X2-X1)**2) ) /Den**2
	Y_Foot= ( (X2-X1)*(Y1*(X2-X0)+Y2*(X0-X1)) + Y0*((Y2-Y1)**2) ) /Den**2


	return Num/Den, X_Foot,Y_Foot




def findDistance(X1,Y1,X2,Y2):
		return np.sqrt((X1-X2)**2+(Y1-Y2)**2)



def getBlockedIntensity(D,Wo):
	p=np.zeros((len(D),1))
	p=0.5*(2-scipy.special.erfc(D/(0.707*Wo)))
	return p


def getPR2BeamSpot(dtheta1DEG):
	global X_spot,Y_spot  
	dtheta1=dtheta1DEG*(np.pi/180)
	X_spot=X_prm+ l1*np.cos(theta1_o+dtheta1)
	Y_spot=Y_prm+ l1*np.sin(theta1_o+dtheta1)

	Diff_X= X_spot-X_pr3	
	Diff_Y= Y_pr3- Y_spot	

	theta2=np.arctan(Diff_Y/Diff_X)

	dtheta2=theta2- theta2_o	

	# print(f'PRM {X_prm:0.02f}, {Y_prm:0.02f}')
	# print(f'PR3 {X_pr3:0.02f}, {Y_pr3:0.02f}')
	# print(f'SPT {X_spot:0.02f}, {Y_spot:0.02f}')

	# print('theta1_o',theta1_o*(180/np.pi))
	# print('theta2_o',theta2_o*(180/np.pi))

	# print('dtheta1',dtheta1*(180/np.pi))
	# print('dtheta2',dtheta2*(180/np.pi))
	return X_spot,Y_spot,dtheta2*(180/np.pi)




def PlotGeom():
	global X_int_dn	,Y_int_dn, X_int_up	,Y_int_up
	plt.figure(1)
	plt.plot([X_prm,X_spot],[Y_prm,Y_spot],marker='o',color='b')
	plt.plot([X_pr3,X_spot],[Y_pr3,Y_spot],marker='o',color='r')
	# plt.plot([X_topEdge,X_bottomEdge],[Y_topEdge,Y_bottomEdge],marker='o',color='c')
	plt.plot(X_int_up	,Y_int_up,marker='*',color='b')
	plt.plot(X_int_dn	,Y_int_dn,marker='*',color='r')
	
	plt.plot(X_Foot_up	,Y_Foot_up,marker='.',color='k')
	plt.plot(X_Foot_dn	,Y_Foot_dn,marker='.',color='k')
	# plt.grid()

	# plt.plot([X_topEdge,X_Foot_up],[Y_topEdge,Y_Foot_up],linestyle='--',color='k')


	# X_Foot,Y_Foot

	# print('Baffle Center', (X_topEdge+X_bottomEdge)/2 ,(Y_topEdge+Y_bottomEdge)/2)
	# plt.gca().set_aspect('equal')





def ScanTheta():
	global X_int_dn	,Y_int_dn, X_int_up	,Y_int_up,Theta,Theta2,X,Y,D1,D2,ND1,ND2,X_Foot_up,Y_Foot_up,X_Foot_dn,Y_Foot_dn

	n=20


	Theta_start=-.105
	Theta_stop=.085
	# Theta_start=-.04
	# Theta_stop=.04

	Theta= np.linspace(Theta_start,Theta_stop,2*n+1)
	
	# Thetamax=.095
	# Theta= Thetamax*np.linspace(-n,n,2*n+1)/n
	

	# print(Theta)
	# print(len(Theta))

	X=np.zeros((2*n+1,))
	Y=np.zeros((2*n+1,))
	Theta2=np.zeros((2*n+1,))
	D1=np.zeros((2*n+1,))
	D2=np.zeros((2*n+1,))

	ND1=np.zeros((2*n+1,))
	ND2=np.zeros((2*n+1,))


	for i in range(2*n+1):
		setGeom(Theta[i])
		X_spot,Y_spot,Theta2[i] = getPR2BeamSpot(Theta[i])
		X_int_up,Y_int_up=findIntersection(X_topEdge,Y_topEdge,X_bottomEdge,Y_bottomEdge,X_spot,Y_spot,X_prm,Y_prm)
		X_int_dn,Y_int_dn=findIntersection(X_topEdge,Y_topEdge,X_bottomEdge,Y_bottomEdge,X_spot,Y_spot,X_pr3,Y_pr3)



		D1[i]=findDistance(X_int_dn,Y_int_dn,X_bottomEdge,Y_bottomEdge)
		D2[i]=findDistance(X_int_up,Y_int_up,X_topEdge,Y_topEdge)
		ND2[i],X_Foot_up,Y_Foot_up=findNormalDistance(X_topEdge,Y_topEdge,X_spot,Y_spot,X_pr3,Y_pr3,)	
		ND1[i],X_Foot_dn,Y_Foot_dn=findNormalDistance(X_bottomEdge,Y_bottomEdge,X_prm,Y_prm,X_spot,Y_spot)	

		# print('D1',D1)
		# print('D2',D2)

		

		X[i]=X_spot
		Y[i]=Y_spot
		# print('X_int_dn',X_int_dn)
		if i%8==0: 
			PlotGeom()




# PlotGeom()
# setGemm()
ScanTheta()

Wo1=6.3
Wo2=7.6
Idn=getBlockedIntensity(ND1,Wo1)
Iup=getBlockedIntensity(ND2,Wo2)

# print('X_Foot_up',X_Foot_up)
# print('Y_Foot_up',Y_Foot_up)
plt.figure(1)
offset=10
plt.arrow(X_topEdge,Y_topEdge,0.5*(X_Foot_up-X_topEdge),0.5*(Y_Foot_up-Y_topEdge),width=5,head_length=3) 
plt.arrow(X_bottomEdge,Y_bottomEdge,0.15*(X_Foot_dn-X_bottomEdge),0.15*(Y_Foot_dn-Y_bottomEdge),width=5,head_length=3) 
# plt.text(X_topEdge,Y_topEdge,0.5*(X_Foot_up-X_topEdge),0.5*(Y_Foot_up-Y_topEdge),width=1) 
plt.text(0.5*(X_Foot_up+X_topEdge)+offset, 0.5*(Y_Foot_up+Y_topEdge), 'D2',size=15)
plt.text(0.5*(X_Foot_dn+X_bottomEdge)-4*offset, 0.5*(Y_Foot_dn+Y_bottomEdge)-offset, 'D1',size=15)









# plt.figure(2)

# plt.plot(Theta,Theta2)

# plt.subplot(2,2,1)
plt.figure(figsize=(10,8))

plt.plot(Theta,X,label='Axial')
plt.plot(Theta,Y,label='Horizontal')
plt.legend()
plt.xlabel('Theta X (Radian)')
plt.ylabel('PR2_Spot_Location')

plt.savefig('Plot1.png')
plt.figure(figsize=(10,8))

# plt.subplot(2,2,3)

# plt.plot(Y,D1,label='D1')
# plt.plot(Y,D2,label='D2')

plt.plot(Y,ND1,label='Lower Edge (D1)')
plt.plot(Y,ND2,label='Upper Edge (D2)')
plt.xlabel('PR2_Spot_Location_Horizontal (mm)')
plt.ylabel('Beam to Baffle Edge Distance (mm)')
plt.legend()
plt.savefig('Plot2.png')


# plt.subplot(2,2,2)
plt.figure(figsize=(10,8))

plt.plot(Y, Idn,label='Lower Edge Transmission')
plt.plot(Y, Iup,label='Upper Edge Transmission')

# plt.semilogy(Y, Idn,label='Idn')
# plt.semilogy(Y, Iup,label='Iup')


plt.xlabel('PR2_Spot_Location_Horizontal (mm)')
plt.ylabel('Transmission')
plt.legend()
plt.savefig('Plot3.png')


# plt.subplot(2,2,4)
plt.figure(figsize=(10,8))

Placement_offset=0
plt.plot(Pr2_Spot_19Feb-Placement_offset, Power_19Feb /np.max(0.00464),label='Expt Data (19 Feb 2025)',alpha=.6)
plt.plot(Expt_PR2_Spots, Expt_CIRC_Powers/np.max(Expt_CIRC_Powers),label='Expt Data (5 July 2024)')

plt.plot(Y, Iup*Idn,label='Simulated Transmission')
# plt.semilogy(Y, Iup*Idn,label='Net Power')
plt.xlabel('PR2_Spot_Location_Horizontal (mm)')

plt.ylabel('Net Transmission')
plt.legend()
plt.savefig('Plot4.png')


plt.show()