!***********************************************
!*  Makro zur Drehung des Elementkoordinaten-  *
!*  systems um die z-Achse fuer Shell181       *
!*                                             *
!*  arg1: Rotationswinkel                      *
!*                                             *
!***********************************************

*afun,deg							!Winkel in Grad
/psymb,cs,0							!lokale KSYS ausblenden
EPS = 1e-5

*get,ACTI_CSYS,active,,csys			!aktives Koordinatensystem
*get,ACTI_ESYS,active,,esys			!aktives Elementkoordinatensystem

cm,ELEM_ROTESYS,elem
*get,E_COUNT,elem,,count			!Anzahl selektierter Elemente
*dim,EARRAY,array,E_COUNT,1
*vget,EARRAY(1,1),elem,,elist		!Elementnummern in Array abspeichern

*do,ii,1,E_COUNT,1

	E_NUM = EARRAY(ii,1)
	
	! Normale (->Dickenrichtung) bestimmen
	!**************************************
	csys,0
	wpcsys,,0
	
	*get,NODE1,elem,E_NUM,node,1
	*get,NODE2,elem,E_NUM,node,2
	*get,NODE3,elem,E_NUM,node,3
	!Knotenkoordinaten auslesen
	a1 = nx(NODE2)-nx(NODE1)
	a2 = ny(NODE2)-ny(NODE1)
	a3 = nz(NODE2)-nz(NODE1)
	b1 = nx(NODE3)-nx(NODE1)
	b2 = ny(NODE3)-ny(NODE1)
	b3 = nz(NODE3)-nz(NODE1)
	!Vektorprodukt bilden
	n1 = a2*b3-a3*b2
	n2 = a3*b1-a1*b3
	n3 = a1*b2-a2*b1
		
	!WorkingPlane ins Element setzen
	esel,s,,,E_NUM
	nsle,s
	nwpave,all
	
	! z-Achse ausrichten
	!********************
	!Drehwinkel um y-Achse
	*if,abs(n3),lt,EPS,then
		PHI_Y = 90
	*else
		PHI_Y = atan(n1/n3)
	*endif
	
	!Drehwinkel um x-Achse
	PHI_X = atan(-n2/(sin(PHI_Y)*n1+cos(PHI_Y)*n3))
	
	!WorkingPlane drehen, sodass x- und y- Achse senkrecht auf Normale
	wprota,,,PHI_Y
	wprota,,PHI_X
	
	! Winkel zur Ausrichtung der x-Achse bestimmen
	!**********************************************
	*get,ELEM_SYS,elem,E_NUM,attr,esys		!Elementkoordinatensystem
	
	*if,ELEM_SYS,eq,0,then	
		zaehler = sin(PHI_X)*sin(PHI_Y)*a1 + cos(PHI_X)*a2 + sin(PHI_X)*cos(PHI_Y)*a3
		nenner = cos(PHI_Y)*a1-sin(PHI_Y)*a3
		PHI_Z = atan(zaehler / nenner)
		*if,abs(a1),gt,EPS,then
			lambda = a1 / (cos(PHI_Y)*cos(PHI_Z)+sin(PHI_X)*sin(PHI_Y)*sin(PHI_Z))
		*elseif,abs(a2),gt,EPS,then
			lambda = a2 / (cos(PHI_X)*sin(PHI_Z))
		*else
			lambda = a3 / (-sin(PHI_Y)*cos(PHI_Z) + sin(PHI_X)*cos(PHI_Y)*sin(PHI_Z))
		*endif
		
	*else
		*get,ELSYS_TYPE,cdsy,ELEM_SYS,attr,kcs	!Typ des Koordinatensystems
												! 0 = kartesisch
												! 1 = zylindrisch
												
		*get,ESYS_ANGXY,cdsy,ELEM_SYS,ang,xy	!Drehwinkel 
		*get,ESYS_ANGYZ,cdsy,ELEM_SYS,ang,yz
		*get,ESYS_ANGZX,cdsy,ELEM_SYS,ang,zx
		*get,ESYS_LOCX,cdsy,ELEM_SYS,loc,x		!ESYS-Ursprung
		*get,ESYS_LOCY,cdsy,ELEM_SYS,loc,y
		*get,ESYS_LOCZ,cdsy,ELEM_SYS,loc,z
	
		!Drehmatrix
		d11 = cos(ESYS_ANGZX)*cos(ESYS_ANGXY)-sin(ESYS_ANGYZ)*sin(ESYS_ANGZX)*sin(ESYS_ANGXY)
		d12 = -cos(ESYS_ANGYZ)*sin(ESYS_ANGXY)
		d13 = sin(ESYS_ANGZX)*cos(ESYS_ANGXY)+sin(ESYS_ANGYZ)*cos(ESYS_ANGZX)*sin(ESYS_ANGXY)
		d21 = cos(ESYS_ANGZX)*sin(ESYS_ANGXY)+sin(ESYS_ANGYZ)*sin(ESYS_ANGZX)*cos(ESYS_ANGXY)
		d22 = cos(ESYS_ANGYZ)*cos(ESYS_ANGXY)
		d23 = sin(ESYS_ANGZX)*sin(ESYS_ANGXY)-sin(ESYS_ANGYZ)*cos(ESYS_ANGZX)*cos(ESYS_ANGXY)
		d31 = -cos(ESYS_ANGYZ)*sin(ESYS_ANGZX)
		d32 = sin(ESYS_ANGYZ)
		d33 = cos(ESYS_ANGYZ)*cos(ESYS_ANGZX)
												
		*if,ELSYS_TYPE,eq,0,then				!kart. ESYS		
			x1 = 1								!x-Achse in kart. ESYS
			x2 = 0
			x3 = 0
			
			y1 = 0								!y-Achse in kart. ESYS
			y2 = 1
			y3 = 0
			
		*elseif,ELSYS_TYPE,eq,1,then			!zyl. ESYS
			x1 = centrx(E_NUM)-ESYS_LOCX		!x-Achse in zyl. ESYS
			x2 = centry(E_NUM)-ESYS_LOCY
			x3 = 0
			
			y1 = -x2							!y-Achse in zyl. ESYS
			y2 = x1
			y3 = 0
			
		*endif
		
		!ESYS-x-Achsen-Vektor im globalen KSYS bestimmen
		x1_glob = d11*x1 + d12*x2 + d13*x3
		x2_glob = d21*x1 + d22*x2 + d23*x3
		x3_glob = d31*x1 + d32*x2 + d33*x3
		
		!ESYS-y-Achsen-Vektor im globalen KSYS
		y1_glob = d11*y1 + d12*y2 + d13*y3
		y2_glob = d21*y1 + d22*y2 + d23*y3
		y3_glob = d31*y1 + d32*y2 + d33*y3		
		
		!Winkel zwischen x-Achse (global) und Normale ermitteln
		zaehler = x1_glob*n1 + x2_glob*n2 + x3_glob*n3
		nenner = sqrt(x1_glob**2 + x2_glob**2 + x3_glob**2)*sqrt(n1**2 + n2**2 + n3**2)
		alpha = acos(zaehler / nenner)
		
		!WP-x-Achse nach ESYS-y
		*if,alpha,le,45,then
			zaehler = sin(PHI_X)*sin(PHI_Y)*y1_glob + cos(PHI_X)*y2_glob + sin(PHI_X)*cos(PHI_Y)*y3_glob
			nenner = cos(PHI_Y)*y1_glob - sin(PHI_Y)*y3_glob
			*if,abs(nenner),lt,EPS,then
				PHI_Z = 90
			*else
				PHI_Z = atan(zaehler / nenner)
			*endif
			
			*if,abs(y1_glob),gt,EPS,then
				lambda = y1_glob / (cos(PHI_Y)*cos(PHI_Z)+sin(PHI_X)*sin(PHI_Y)*sin(PHI_Z))
			*elseif,abs(y2_glob),gt,EPS,then
				lambda = y2_glob / (cos(PHI_X)*sin(PHI_Z))
			*else
				lambda = y3_glob / (-sin(PHI_Y)*cos(PHI_Z) + sin(PHI_X)*cos(PHI_Y)*sin(PHI_Z))
			*endif
		
		!WP-x-Achse nach ESYS-x
		*else
			zaehler = sin(PHI_X)*sin(PHI_Y)*x1_glob + cos(PHI_X)*x2_glob + sin(PHI_X)*cos(PHI_Y)*x3_glob
			nenner = cos(PHI_Y)*x1_glob - sin(PHI_Y)*x3_glob
			*if,abs(nenner),lt,EPS,then
				PHI_Z = 90
			*else
				PHI_Z = atan(zaehler / nenner)
			*endif

			*if,abs(x1_glob),gt,EPS,then
				lambda = x1_glob / (cos(PHI_Y)*cos(PHI_Z)+sin(PHI_X)*sin(PHI_Y)*sin(PHI_Z))
			*elseif,abs(x2_glob),gt,EPS,then
				lambda = x2_glob / (cos(PHI_X)*sin(PHI_Z))
			*else
				lambda = x3_glob / (-sin(PHI_Y)*cos(PHI_Z) + sin(PHI_X)*cos(PHI_Y)*sin(PHI_Z))
			*endif
			
		*endif
	*endif
	
	! x-Achse durch Drehung um z ausrichten, ggf. Richtungssinn anpassen	
	*if,lambda,gt,0,then
		wprota,PHI_Z
	*else
		wprota,PHI_Z + 180
	*endif

	!**************************
	!*  Neues ESYS erstellen  *
	!**************************
	wprota,arg1					!Working Plane um z-Achse drehen
	*get,MAX_CS,cdsy,,num,max
	*if,MAX_CS,lt,10,then
		MAX_CS = 10
	*endif
	cswpla,MAX_CS+1,0			!kart. KSYS in WorkingPlane setzen
	esys,MAX_CS+1
	emodif,E_NUM,esys

*enddo

*del,EARRAY

esel,s,,,ELEM_ROTESYS
nsle,s
cmdele,ELEM_ROTESYS

csys,ACTI_CSYS
wpcsys,,ACTI_CSYS
esys,ACTI_ESYS

*afun,rad