! File: submodule_lorentz_group_constructors.f90 ! Authors: Francesco Torsello (FT) !************************************************************************ ! Copyright (C) 2020-2023 Francesco Torsello * ! * ! This file is part of SPHINCS_ID * ! * ! SPHINCS_ID is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by * ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! SPHINCS_ID is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with SPHINCS_ID. If not, see <https://www.gnu.org/licenses/>. * ! The copy of the GNU General Public License should be in the file * ! 'COPYING'. * !************************************************************************ SUBMODULE(lorentz_group) constructors !*********************************************************** ! !# This submodule contains the implementations of the ! constructors of the TYPES representing the members of ! the Lorentz group defined in MODULE lorentz_group, ! presently boosts and spatial rotations ! ! FT 08.12.2022 ! !*********************************************************** USE utility, ONLY: zero, one IMPLICIT NONE CONTAINS !------------! !-- BOOSTS --! !------------! MODULE PROCEDURE construct_boost !*********************************************************** ! !# Construct a boost transformation from the spatial velocity ! ! FT 08.12.2022 ! !*********************************************************** IMPLICIT NONE INTEGER:: i, j DOUBLE PRECISION:: v_sqnorm, delta DOUBLE PRECISION, DIMENSION(4,4):: identity(0:3,0:3) LOGICAL:: is_identity v_sqnorm= euclidean_inner_product(v,v) boost% v_speed= SQRT(v_sqnorm) IF( boost% v_speed >= one )THEN PRINT *, "** ERROR! The speed given to construct_boost is larger", & " than, or equal to, 1!" PRINT *, "|v|^2=", boost% v_speed PRINT *, "* Stopping..." PRINT * STOP ENDIF boost% v= v boost% lambda= one/SQRT( one - v_sqnorm ) boost% p= boost% lambda*boost% v CALL boost% compute_boost_matrices(boost% p, boost% lambda_s, boost% matrix) CALL boost% compute_boost_matrices(- boost% p, boost% inv_lambda_s, & boost% inv_matrix) boost% tr_matrix= TRANSPOSE(boost% matrix) identity= square_matrix_multiplication(boost% inv_matrix, boost% matrix) is_identity= .TRUE. DO i= 0, 3, 1 DO j= 0, 3, 1 IF( i == j )THEN delta= one ELSE delta= zero ENDIF is_identity= is_identity .AND. (identity(i,j) - delta) < 1.D-12 ENDDO ENDDO IF( .NOT.is_identity )THEN PRINT *, "** ERROR! boost% inv_matrix is not the inverse of ", & "boost% matrix! Their product is:" PRINT *, " id(0,:)=", identity(0,:) PRINT *, " id(1,:)=", identity(1,:) PRINT *, " id(2,:)=", identity(2,:) PRINT *, " id(3,:)=", identity(3,:) PRINT * STOP ENDIF END PROCEDURE construct_boost MODULE PROCEDURE construct_boost_components !*********************************************************** ! !# Construct a boost transformation from the components ! of the spatial velocity ! ! FT 10.12.2022 ! !*********************************************************** IMPLICIT NONE boost= lorentz_boost([vx,vy,vz]) END PROCEDURE construct_boost_components MODULE PROCEDURE compute_boost_matrices !*********************************************************** ! !# Computes the spatial part of the matrix of the Lorentz ! boost, and its whole matrix, starting from the vector ! \(p\) ! ! FT 09.12.2022 ! !*********************************************************** USE utility, ONLY: one IMPLICIT NONE DOUBLE PRECISION:: lambda, lambda_plus_one lambda= SQRT(one + euclidean_inner_product(p,p)) lambda_plus_one= lambda + one lambda_s(1)= one + (p(1)**2) /lambda_plus_one lambda_s(2)= (p(1)*p(2))/lambda_plus_one lambda_s(3)= (p(1)*p(3))/lambda_plus_one lambda_s(4)= one + (p(2)**2) /lambda_plus_one lambda_s(5)= (p(2)*p(3))/lambda_plus_one lambda_s(6)= one + (p(3)**2) /lambda_plus_one matrix(0,0) = lambda matrix(0,1:3)= p matrix(1:3,0)= p matrix(1,1:3)= lambda_s(1:3) matrix(2,1:3)= [lambda_s(2), lambda_s(4), lambda_s(5)] matrix(3,1:3)= [lambda_s(3), lambda_s(5), lambda_s(6)] END PROCEDURE compute_boost_matrices !------------------------! !-- SPATIAL ROTATIONS --! !------------------------! MODULE PROCEDURE construct_rotation !*********************************************************** ! !# Construct a spatial rotation from a vector containing ! the Euler angles ! ! FT 09.12.2022 ! !*********************************************************** IMPLICIT NONE INTEGER:: i, j DOUBLE PRECISION:: v_sqnorm, delta DOUBLE PRECISION, DIMENSION(4,4):: identity(0:3,0:3) LOGICAL:: is_identity rotation% euler_angles= euler_angles CALL rotation% compute_rotation_matrices(rotation% euler_angles, & rotation% r_x, rotation% r_y, rotation% r_z, & rotation% r, rotation% matrix, & rotation% inv_r_x, rotation% inv_r_y, rotation% inv_r_z, & rotation% inv_r, rotation% inv_matrix) rotation% tr_r= TRANSPOSE(rotation% r) rotation% tr_matrix= TRANSPOSE(rotation% matrix) identity= square_matrix_multiplication(rotation% inv_matrix, & rotation% matrix) is_identity= .TRUE. DO i= 0, 3, 1 DO j= 0, 3, 1 IF( i == j )THEN delta= one ELSE delta= zero ENDIF is_identity= is_identity .AND. (identity(i,j) - delta) < 1.D-12 ENDDO ENDDO IF( .NOT.is_identity )THEN PRINT *, "** ERROR! rotation% inv_matrix is not the inverse of ", & "rotation% matrix! Their product is:" PRINT *, " id(0,:)=", identity(0,:) PRINT *, " id(1,:)=", identity(1,:) PRINT *, " id(2,:)=", identity(2,:) PRINT *, " id(3,:)=", identity(3,:) PRINT * STOP ENDIF is_identity= .TRUE. DO i= 1, 3, 1 DO j= 1, 3, 1 is_identity= is_identity .AND. & (rotation% inv_r(i,j) - rotation% tr_r(i,j)) < 1.D-12 ENDDO ENDDO IF( .NOT.is_identity )THEN PRINT *, "** ERROR! rotation% inv_r is not the same as ", & "rotation% tr_r! " PRINT * STOP ENDIF END PROCEDURE construct_rotation MODULE PROCEDURE construct_rotation_angles !*********************************************************** ! !# Construct a spatial rotation from the Euler angles ! ! FT 10.12.2022 ! !*********************************************************** IMPLICIT NONE rotation= spatial_rotation([alpha,beta,gamma]) END PROCEDURE construct_rotation_angles MODULE PROCEDURE compute_rotation_matrices !*********************************************************** ! !# Computes the spatial part of the matrix of the Lorentz ! rotation, and its whole matrix, starting from the vector ! \(p\) ! ! FT 09.12.2022 ! !*********************************************************** IMPLICIT NONE r_x(1,:)= [one, zero, zero] r_x(2,:)= [zero, COS(euler_angles(1)), - SIN(euler_angles(1))] r_x(3,:)= [zero, SIN(euler_angles(1)), COS(euler_angles(1))] r_y(1,:)= [ COS(euler_angles(2)), zero, SIN(euler_angles(2))] r_y(2,:)= [ zero, one, zero] r_y(3,:)= [- SIN(euler_angles(2)), zero, COS(euler_angles(2))] r_z(1,:)= [COS(euler_angles(3)), - SIN(euler_angles(3)), zero] r_z(2,:)= [SIN(euler_angles(3)), COS(euler_angles(3)), zero] r_z(3,:)= [zero, zero, one] r= square_matrix_multiplication(r_z,square_matrix_multiplication(r_y,r_x)) matrix(0,0)= one matrix(0,1:3)= [zero,zero,zero] matrix(1:3,0)= [zero,zero,zero] matrix(1,1:3)= r(1,1:3) matrix(2,1:3)= r(2,1:3) matrix(3,1:3)= r(3,1:3) ! !-- Computation of the inverse matrices ! inv_r_x(1,:)= [one, zero, zero] inv_r_x(2,:)= [zero, COS(-euler_angles(1)), - SIN(-euler_angles(1))] inv_r_x(3,:)= [zero, SIN(-euler_angles(1)), COS(-euler_angles(1))] inv_r_y(1,:)= [ COS(-euler_angles(2)), zero, SIN(-euler_angles(2))] inv_r_y(2,:)= [ zero, one, zero] inv_r_y(3,:)= [- SIN(-euler_angles(2)), zero, COS(-euler_angles(2))] inv_r_z(1,:)= [COS(-euler_angles(3)), - SIN(-euler_angles(3)), zero] inv_r_z(2,:)= [SIN(-euler_angles(3)), COS(-euler_angles(3)), zero] inv_r_z(3,:)= [zero, zero, one] inv_r= square_matrix_multiplication(inv_r_x, & square_matrix_multiplication(inv_r_y,inv_r_z)) inv_matrix(0,0)= one inv_matrix(0,1:3)= [zero,zero,zero] inv_matrix(1:3,0)= [zero,zero,zero] inv_matrix(1,1:3)= inv_r(1,1:3) inv_matrix(2,1:3)= inv_r(2,1:3) inv_matrix(3,1:3)= inv_r(3,1:3) END PROCEDURE compute_rotation_matrices MODULE PROCEDURE get_lambda !*********************************************************** ! !# Returns the Lorentz factor [[lorentz_boost:lambda]] ! ! FT 21.02.2023 ! !*********************************************************** IMPLICIT NONE lambda= this% lambda END PROCEDURE get_lambda END SUBMODULE constructors