! lcsx=x3a ! lcsz=VecNormalize(x3av12a) ! lcsy=lcszx3a
!
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(vector), | dimension(:) | :: | Surf | |||
| integer, | intent(in) | :: | NSides | |||
| real(kind=r64), | intent(out) | :: | Azimuth | |||
| real(kind=r64), | intent(out) | :: | Tilt | |||
| type(vector) | :: | lcsx | ||||
| type(vector) | :: | lcsy | ||||
| type(vector) | :: | lcsz | ||||
| real(kind=r64), | intent(in) | :: | surfaceArea | |||
| type(vector), | intent(in) | :: | NewellSurfaceNormalVector |
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Nodes of different colours represent the following:
Solid arrows point from a procedure to one which it calls. Dashed arrows point from an interface to procedures which implement that interface. This could include the module procedures in a generic interface or the implementation in a submodule of an interface in a parent module. Where possible, edges connecting nodes are given different colours to make them easier to distinguish in large graphs.
Subroutine DetermineAzimuthAndTilt(Surf,NSides,Azimuth,Tilt,lcsx,lcsy,lcsz,surfaceArea,NewellSurfaceNormalVector)
! PURPOSE OF THIS SUBROUTINE:
! This subroutine determines the Azimuth (outward normal) angle,
! Tilt angle of a given surface defined by the set of input vectors.
! REFERENCE:
! Discussions and examples from Bill Carroll, LBNL.
IMPLICIT NONE
! SUBROUTINE ARGUMENT DEFINITIONS:
type (vector), dimension(:) :: Surf ! Surface Definition
integer, intent(in) :: NSides ! Number of sides to surface
real(r64), intent(out) :: Azimuth ! Outward Normal Azimuth Angle
real(r64), intent(out) :: Tilt ! Tilt angle of surface
type (vector) :: lcsx
type (vector) :: lcsy
type (vector) :: lcsz
real(r64), intent(in) :: surfaceArea
type (vector), intent(in) :: NewellSurfaceNormalVector
character(len=*), parameter :: fmt3="(A,3(1x,f18.13))"
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
! type(vector) :: x3,y3,z3,v12
type(vector) :: x2
type(vector) :: x3a,v12a
type(vector) :: v0
type(vector) :: v1
type(vector) :: v2
type(vector) :: cs3_2
type(vector) :: cs3_0
type(vector) :: cs3_1
! type(vector) :: y2
type(vector) :: z3
real(r64) costheta
real(r64) rotang_0
! real(r64) rotang_2
real(r64) az
! real(r64) azm
real(r64) tlt
! real(r64) newtlt
! real(r64) roundval
! real(r64) xcomp
! real(r64) ycomp
! real(r64) zcomp
! real(r64) proj
! integer :: scount
! integer :: nvert1
! real(r64) :: tltcos
!!! x3=VecNormalize(Surf(2)-Surf(1))
!!! v12=Surf(3)-Surf(2)
!!!
!!! z3=VecNormalize(x3*v12)
!!! y3=z3*x3
!!! roundval=10000.d0
!!! CALL VecRound(x3,roundval)
!!! CALL VecRound(y3,roundval)
!!! CALL VecRound(z3,roundval)
!!!
!!!! Direction cosines, local coordinates.
!!!! write(OUTPUT,*) 'lcs:'
!!!! write(OUTPUT,*) 'x=',x3
!!!! write(OUTPUT,*) 'y=',y3
!!!! write(OUTPUT,*) 'z=',z3
x3a=VecNormalize(Surf(3)-Surf(2))
v12a=Surf(1)-Surf(2)
!!! lcsx=x3a
!!! lcsz=VecNormalize(x3a*v12a)
!!! lcsy=lcsz*x3a
lcsx=x3a
lcsz=NewellSurfaceNormalVector
lcsy=lcsz*x3a
!!!
! Vec3d v0(p1 - p0); ! BGL has different conventions...p0=surf(2), etc
v0=Surf(3)-Surf(2)
! Vec3d v1(p2 - p0);
v1=Surf(1)-Surf(2)
! Vec3d v2 = cross(v0,v1);
v2=v0*v1
! cs3[2] = norm(v2); // z
cs3_2=VecNormalize(v2)
! cs3[0] = norm(v0); // x
cs3_0=VecNormalize(v0)
! cs3[1] = cross(cs3[2],cs3[0]); // y
cs3_1=cs3_2*cs3_0
! Vec3d z3 = cs3[2];
z3=cs3_2
! double costheta = dot(z3,Ref_CS[2]);
costheta = z3 .dot. zunit
! if ( fabs(costheta) < 1.0d0) { // normal cases
if (abs(costheta) < 1.0d0) then
! // azimuth
! Vec3d x2 = cross(Ref_CS[2],z3); // order is important; x2 = x1
! RotAng[0] = atan2(dot(x2,Ref_CS[1]),dot(x2,Ref_CS[0]));
x2= zunit*z3
rotang_0 = atan2(x2 .dot. yunit, x2 .dot. xunit)
else
! }
! else { // special cases: tilt angle theta = 0, PI
! // azimuth
! RotAng[0] = atan2(dot(cs3[0],Ref_CS[1]),dot(cs3[0],Ref_CS[0]) );
rotang_0=atan2(cs3_0 .dot. yunit, cs3_0 .dot. xunit)
!
endif
tlt=acos(NewellSurfaceNormalVector%z)
tlt=tlt/degtoradians
az=rotang_0
az=az/degtoradians
az=MOD(450.d0 - az, 360.d0)
az=az+90.d0
if (az < 0.0d0) az=az+360.d0
az=mod(az,360.d0)
! Normalize the azimuth angle so it is positive
if (abs(az-360.d0) < 1.d-3) az=0.0d0
Azimuth=az
Tilt=tlt
RETURN
END Subroutine DetermineAzimuthAndTilt