I have a coordinate system XYZ (black on figure) which can be rotated around
Z axis (Zrot angle) or around Y axis (Yrot angle) or around both axisies (around Z first and then around rotated (in first rotation around Z axis) Y axis).
E.g. after rotation only around Z axis we get X'Y'Z coordinate system (see fig).
I have also the direction (point P) defined by azimuth and evation (e.g. AZ, EL) in XYZ coordinate system (blue on figure).
Rotation code are below:
% geodetic elevation and azimuth in deg(El, AZ)
Esvo = 50;
Asvo = 30;
% MatLab elevation and azimuth 
Esv = Esvo;
Asv = Asvo;
 if Asvo < 180; Asv=-Asvo;
 else Asv = 360-Asvo;
 end
r=1
dtr = pi/180;
Es=Esv*dtr  % deg to rad conversion
As=Asv*dtr  % deg to rad conversion
% initial position in MatLab
[x0,y0,z0] = sph2cart(As,Es,r)
poz0 = [x0 y0 z0]
% coordinate system rotation angles in deg (Zrot=Aro1; Yrot=Ero1)
    Ero1 = 10;
    Aro1 = 40;
    
% AZ change after Zrot rotation (geodetic)
	Asvn = Asvo-Aro1	
	% chenged AZ in MatLab
	 if Asvn < 180; Asvn=-Asvn;
	 else Asvn = 360-Asvn;
	 end
Esv=Esv*dtr  % deg to rad conversion
Asvn=Asvn*dtr  % deg to rad conversion
% position after rotation around Z-axis	 
	[x1,y1,z1] = sph2cart(Asvn,Esv,r)
	poz1 = [x1 y1 z1]
% rotation around Y-axis 
    ROTY = [ cosd(Ero1) 0  sind(Ero1);
                  0     1        0;
            -sind(Ero1) 0  cosd(Ero1)];
  
    % position after two rotation (in MatLab)
    poz2 = poz1*ROTY;
    x2 = poz2(1);
    y2 = poz2(2);
    z2 = poz2(3);
	
I would like to create a PLANE (red on figure) which is perpendicular to direction P (direction P will be normal to PLANE) and then calculate angle B
which is the difference betwen X and X' axisies orthogonally projected on PLANE (red on fig).
Any suggestions how can be it done?