From 1278c94c2c5903b9543f1c1d0a19a7bb02126dde Mon Sep 17 00:00:00 2001 From: Behdad Shaarbaf Ebrahimi Date: Wed, 3 Dec 2025 16:11:33 +1300 Subject: [PATCH 1/3] export file completed will WSS transient values --- src/lib/wave_transmission.f90 | 1131 +++++++++++++++++---------------- 1 file changed, 567 insertions(+), 564 deletions(-) diff --git a/src/lib/wave_transmission.f90 b/src/lib/wave_transmission.f90 index 40332ae4..e4b536c4 100644 --- a/src/lib/wave_transmission.f90 +++ b/src/lib/wave_transmission.f90 @@ -74,11 +74,11 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, real(dp), allocatable :: reflected_pressure_previous(:) real(dp), allocatable :: forward_flow(:) real(dp), allocatable :: reflected_flow(:) - real(dp), allocatable :: terminals_radius(:),WSS(:) - real(dp), allocatable :: p_terminal(:),p_previous(:),terminal_flow(:) + real(dp), allocatable :: terminals_radius(:),WSS(:),elem_rad(:),alpha(:) + real(dp), allocatable :: p_terminal(:),p_previous(:),terminal_flow(:),elem_flow(:) integer :: min_art,max_art,min_ven,max_ven,min_cap,max_cap,ne,nu,nt,nf,np,np_previous,ne_previous character(len=30) :: tree_direction,mechanics_type - real(dp) start_time,end_time,dt,time,omega,delta_p + real(dp) start_time,end_time,dt,time,omega,delta_p,Ppl real(dp) grav_vect(3),grav_factor,mechanics_parameters(2) integer :: AllocateStatus,fid=10,fid2=20,fid3=30,fid4=40,fid5=50,fid6=60,fid7=70,fid8=80,fid9=90 integer :: fid10=100,fid11=110,fid12=120,fid13=130 @@ -237,12 +237,19 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, if (AllocateStatus /= 0) STOP "*** Not enough memory for terminals_radius array ***" allocate (WSS(n_time)) if (AllocateStatus /= 0) STOP "*** Not enough memory for WSS array ***" + allocate (elem_rad(n_time)) + if (AllocateStatus /= 0) STOP "*** Not enough memory for elem_rad array ***" + allocate (elem_flow(n_time)) + if (AllocateStatus /= 0) STOP "*** Not enough memory for elem_flow array ***" + allocate (alpha(num_elems)) + if (AllocateStatus /= 0) STOP "*** Not enough memory for elasticity array ***" !initialise admittance char_admit=0.0_dp eff_admit=0.0_dp + alpha=0.0_dp !calculate characteristic admittance of each branch - call characteristic_admittance(no_freq,char_admit,prop_const,harmonic_scale, & + call characteristic_admittance(no_freq,char_admit,prop_const,alpha,harmonic_scale, & density,viscosity,admit_param,elast_param,mechanics_parameters,grav_vect,remodeling_grade) !Apply boundary conditions to terminal units @@ -720,203 +727,203 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, write(10, *) " ""LBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," write(10, *) " ""unit"": ""radians""" write(10, *) " }," - write(10, *) " ""Flow phase"":{" - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Calculating LUL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,11)),real(q_factor(i,11), 8))+b(i) - enddo - write(10, *) " ""LUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+10)),real(q_factor(i,min_ven+10), 8))+b(i) - enddo - write(10, *) " ""LUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,20)),real(q_factor(i,20), 8))+b(i) - enddo - write(10, *) " ""LLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+19)),real(q_factor(i,min_ven+19), 8))+b(i) - enddo - write(10, *) " ""LLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RUL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,15)),real(q_factor(i,15), 8))+b(i) - enddo - write(10, *) " ""RUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+14)),real(q_factor(i,min_ven+14), 8))+b(i) - enddo - write(10, *) " ""RUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,23)),real(q_factor(i,23), 8))+b(i) - enddo - write(10, *) " ""RLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+22)),real(q_factor(i,min_ven+22), 8))+b(i) - enddo - write(10, *) " ""RLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RML_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,24)),real(q_factor(i,24), 8))+b(i) - enddo - write(10, *) " ""RML_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RML_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+23)),real(q_factor(i,min_ven+23), 8))+b(i) - enddo - write(10, *) " ""RML_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,1)),real(q_factor(i,1), 8))+b(i)!-atan2(dimag(char_admit(i,1))& - !,real(char_admit(i,1), 8)) - enddo - write(10, *) " ""MPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,2)),real(q_factor(i,2), 8))+b(i) - enddo - write(10, *) " ""MPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,3)),real(q_factor(i,3), 8))+b(i) - enddo - write(10, *) " ""MPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,4)),real(q_factor(i,4), 8))+b(i) - enddo - write(10, *) " ""MPA_A_4"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,5)),real(q_factor(i,5), 8))+b(i) - enddo - write(10, *) " ""LPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,6)),real(q_factor(i,6), 8))+b(i) - enddo - write(10, *) " ""LPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,7)),real(q_factor(i,7), 8))+b(i) - enddo - write(10, *) " ""LPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,8)),real(q_factor(i,8), 8))+b(i) - enddo - write(10, *) " ""RPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,9)),real(q_factor(i,9), 8))+b(i) - enddo - write(10, *) " ""RPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,10)),real(q_factor(i,10), 8))+b(i) - enddo - write(10, *) " ""RPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,16)),real(q_factor(i,16), 8))+b(i) - enddo - write(10, *) " ""RBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+15)),real(q_factor(i,min_ven+15), 8))+b(i) - enddo - write(10, *) " ""RBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,12)),real(q_factor(i,12), 8))+b(i) - enddo - write(10, *) " ""LBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = 0 - do i = 1, no_freq - imped(i+1) = atan2(dimag(q_factor(i,min_ven+11)),real(q_factor(i,min_ven+11), 8))+b(i) - enddo - write(10, *) " ""LBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - write(10, *) " ""unit"": ""radians""" - write(10, *) " }," + ! write(10, *) " ""Flow phase"":{" + ! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! ! Calculating LUL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,11)),real(q_factor(i,11), 8))+b(i) + ! enddo + ! write(10, *) " ""LUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+10)),real(q_factor(i,min_ven+10), 8))+b(i) + ! enddo + ! write(10, *) " ""LUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,20)),real(q_factor(i,20), 8))+b(i) + ! enddo + ! write(10, *) " ""LLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+19)),real(q_factor(i,min_ven+19), 8))+b(i) + ! enddo + ! write(10, *) " ""LLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RUL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,15)),real(q_factor(i,15), 8))+b(i) + ! enddo + ! write(10, *) " ""RUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+14)),real(q_factor(i,min_ven+14), 8))+b(i) + ! enddo + ! write(10, *) " ""RUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,23)),real(q_factor(i,23), 8))+b(i) + ! enddo + ! write(10, *) " ""RLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+22)),real(q_factor(i,min_ven+22), 8))+b(i) + ! enddo + ! write(10, *) " ""RLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RML_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,24)),real(q_factor(i,24), 8))+b(i) + ! enddo + ! write(10, *) " ""RML_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RML_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+23)),real(q_factor(i,min_ven+23), 8))+b(i) + ! enddo + ! write(10, *) " ""RML_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,1)),real(q_factor(i,1), 8))+b(i)!-atan2(dimag(char_admit(i,1))& + ! !,real(char_admit(i,1), 8)) + ! enddo + ! write(10, *) " ""MPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,2)),real(q_factor(i,2), 8))+b(i) + ! enddo + ! write(10, *) " ""MPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,3)),real(q_factor(i,3), 8))+b(i) + ! enddo + ! write(10, *) " ""MPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,4)),real(q_factor(i,4), 8))+b(i) + ! enddo + ! write(10, *) " ""MPA_A_4"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,5)),real(q_factor(i,5), 8))+b(i) + ! enddo + ! write(10, *) " ""LPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,6)),real(q_factor(i,6), 8))+b(i) + ! enddo + ! write(10, *) " ""LPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,7)),real(q_factor(i,7), 8))+b(i) + ! enddo + ! write(10, *) " ""LPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,8)),real(q_factor(i,8), 8))+b(i) + ! enddo + ! write(10, *) " ""RPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,9)),real(q_factor(i,9), 8))+b(i) + ! enddo + ! write(10, *) " ""RPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,10)),real(q_factor(i,10), 8))+b(i) + ! enddo + ! write(10, *) " ""RPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,16)),real(q_factor(i,16), 8))+b(i) + ! enddo + ! write(10, *) " ""RBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+15)),real(q_factor(i,min_ven+15), 8))+b(i) + ! enddo + ! write(10, *) " ""RBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,12)),real(q_factor(i,12), 8))+b(i) + ! enddo + ! write(10, *) " ""LBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = 0 + ! do i = 1, no_freq + ! imped(i+1) = atan2(dimag(q_factor(i,min_ven+11)),real(q_factor(i,min_ven+11), 8))+b(i) + ! enddo + ! write(10, *) " ""LBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! write(10, *) " ""unit"": ""radians""" + ! write(10, *) " }," write(10, *) " ""Pressure phase"":{" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1113,202 +1120,202 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, write(10, *) " ""LBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," write(10, *) " ""unit"": ""radians""" write(10, *) " }," - write(10, *) " ""Flow amplitude"":{" - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - ! Calculating LUL_A Flow ampl !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,11) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,11))*a(i) - enddo - write(10, *) " ""LUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+10) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+10))*a(i) - enddo - write(10, *) " ""LUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,20) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,20))*a(i) - enddo - write(10, *) " ""LLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+19) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+19))*a(i) - enddo - write(10, *) " ""LLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RUL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,15) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,15))*a(i) - enddo - write(10, *) " ""RUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+14) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+14))*a(i) - enddo - write(10, *) " ""RUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,23) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,23))*a(i) - enddo - write(10, *) " ""RLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+22) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+22))*a(i) - enddo - write(10, *) " ""RLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RML_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,24) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,24))*a(i) - enddo - write(10, *) " ""RML_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RML_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+23) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+23))*a(i) - enddo - write(10, *) " ""RML_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,1) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i, 1))*a(i) - enddo - write(10, *) " ""MPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,2) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i, 2))*a(i) - enddo - write(10, *) " ""MPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,3) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i, 3))*a(i) - enddo - write(10, *) " ""MPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,4) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i, 4))*a(i) - enddo - write(10, *) " ""MPA_A_4"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,5) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,5))*a(i) - enddo - write(10, *) " ""LPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,6) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,6))*a(i) - enddo - write(10, *) " ""LPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,7) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,7))*a(i) - enddo - write(10, *) " ""LPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,8) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,8))*a(i) - enddo - write(10, *) " ""RPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,9) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,9))*a(i) - enddo - write(10, *) " ""RPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,10) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,10))*a(i) - enddo - write(10, *) " ""RPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,16) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,16))*a(i) - enddo - write(10, *) " ""RBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating RBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+15) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+15))*a(i) - enddo - write(10, *) " ""RBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,12) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,12))*a(i) - enddo - write(10, *) " ""LBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !Calculating LBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - imped(1) = elem_field(ne_Qdot,min_ven+11) - do i = 1, no_freq - imped(i+1) = abs(q_factor(i,min_ven+11))*a(i) - enddo - write(10, *) " ""LBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," - write(10, *) " ""unit"": ""mm3/s""" - write(10, *) " }," + ! write(10, *) " ""Flow amplitude"":{" + ! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! ! Calculating LUL_A Flow ampl !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,11) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,11))*a(i)!*abs(char_admit(i,11))*a(i) + ! enddo + ! write(10, *) " ""LUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+10) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+10))*a(i)!*abs(char_admit(i,min_ven+10))*a(i) + ! enddo + ! write(10, *) " ""LUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,20) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,20))*a(i)!*abs(char_admit(i,20))*a(i) + ! enddo + ! write(10, *) " ""LLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+19) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+19))*a(i)!*abs(char_admit(i,min_ven+19))*a(i) + ! enddo + ! write(10, *) " ""LLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RUL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,15) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,15))*a(i)!*abs(char_admit(i,15))*a(i) + ! enddo + ! write(10, *) " ""RUL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RUL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+14) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+14))*a(i)!*abs(char_admit(i,min_ven+14))*a(i) + ! enddo + ! write(10, *) " ""RUL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RLL_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,23) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,23))*a(i)!*abs(char_admit(i,23))*a(i) + ! enddo + ! write(10, *) " ""RLL_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RLL_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+22) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+22))*a(i)!*abs(char_admit(i,min_ven+22))*a(i) + ! enddo + ! write(10, *) " ""RLL_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RML_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,24) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,24))*a(i)!*abs(char_admit(i,24))*a(i) + ! enddo + ! write(10, *) " ""RML_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RML_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+23) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+23))*a(i)!*abs(char_admit(i,min_ven+23))*a(i) + ! enddo + ! write(10, *) " ""RML_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,1) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i, 1))*a(i)!*abs(char_admit(i,1))*a(i) + ! enddo + ! write(10, *) " ""MPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,2) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,2))*a(i)!*abs(char_admit(i,2))*a(i) + ! enddo + ! write(10, *) " ""MPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,3) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,3))*a(i)!*abs(char_admit(i,3))*a(i) + ! enddo + ! write(10, *) " ""MPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating MPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,4) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,4))*a(i)!*abs(char_admit(i,4))*a(i) + ! enddo + ! write(10, *) " ""MPA_A_4"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,5) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,5))*a(i)!*abs(char_admit(i,5))*a(i) + ! enddo + ! write(10, *) " ""LPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,6) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,6))*a(i)!*abs(char_admit(i,6))*a(i) + ! enddo + ! write(10, *) " ""LPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,7) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,7))*a(i)!*abs(char_admit(i,7))*a(i) + ! enddo + ! write(10, *) " ""LPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,8) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,8))*a(i)!*abs(char_admit(i,8))*a(i) + ! enddo + ! write(10, *) " ""RPA_A_1"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,9) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,9))*a(i)!*abs(char_admit(i,9))*a(i) + ! enddo + ! write(10, *) " ""RPA_A_2"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RPA phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,10) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,10))*a(i)!*abs(char_admit(i,10))*a(i) + ! enddo + ! write(10, *) " ""RPA_A_3"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,16) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,16))*a(i)!*abs(char_admit(i,16))*a(i) + ! enddo + ! write(10, *) " ""RBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating RBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+15) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+15))*a(i)!*abs(char_admit(i,min_ven+15))*a(i) + ! enddo + ! write(10, *) " ""RBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LBS_A phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,12) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,12))*a(i)!*abs(char_admit(i,12))*a(i) + ! enddo + ! write(10, *) " ""LBS_A"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !Calculating LBS_V phase !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! imped(1) = elem_field(ne_Qdot,min_ven+11) + ! do i = 1, no_freq + ! imped(i+1) = abs(q_factor(i,min_ven+11))*a(i)!*abs(char_admit(i,min_ven+11))*a(i) + ! enddo + ! write(10, *) " ""LBS_V"": [", imped(1), ",", (imped(i),",",i=2,no_freq), imped(no_freq+1), "]," + ! write(10, *) " ""unit"": ""mm3/s""" + ! write(10, *) " }," write(10, *) " ""Pressure amplitude"":{" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1930,10 +1937,13 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, !consider first pressure and flow into the vessel (at x=0) open(fid, file = 'incident_pressure.txt', action='write') open(fid2, file = 'incident_flow.txt',action='write') - open(fid3, file = 'total_pressure.txt',action='write') - open(fid4, file = 'total_flow.txt',action='write') + open(fid3, file = 'terminal_pressure.txt',action='write') + open(fid4, file = 'terminal_flow.txt',action='write') open(fid6, file = 'terminal_radii.txt', action='write') - open(fid7, file = 'WSS.txt',action='write') + open(fid7, file = 'terminal_WSS.txt',action='write') + open(200, file = 'WSS_arterial.txt', action='write') + open(210, file = 'radii.txt', action='write') + open(220, file = 'flows.txt', action='write') do nu =1,num_units ne=units(nu) ! terminal elements ne_previous=elem_cnct(-1,1,ne) @@ -1946,93 +1956,44 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, reflected_flow=0.0_dp terminals_radius=0.0_dp WSS=0.0_dp - if (bc_type.eq.'pressure') then - do nt=1,n_time - do nf=1,no_freq - omega=2*pi*nf*harmonic_scale - forward_pressure(nt)=forward_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*cos(omega*time+b(nf)+& - atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))) - forward_pressure_previous(nt)=forward_pressure_previous(nt)+abs(p_factor(nf,ne_previous))*& - a(nf)*cos(omega*time+b(nf)+atan2(dimag(p_factor(nf,ne_previous)),real(p_factor(nf,ne_previous), 8))) - - reflected_pressure(nt)=reflected_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*& - abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& - cos(omega*time+b(nf)+& - atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& - (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& - atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))) - reflected_pressure_previous(nt)=reflected_pressure_previous(nt)+abs(p_factor(nf,ne_previous))*& - a(nf)*abs(reflect(nf,ne_previous))*exp((-2*elem_field(ne_length,ne_previous))*& - (real(prop_const(nf,ne_previous), 8)))*cos(omega*time+b(nf)+& - atan2(dimag(p_factor(nf,ne_previous)),real(p_factor(nf,ne_previous), 8))+& - (-2*elem_field(ne_length,ne_previous))*(dimag(prop_const(nf,ne_previous)))+& - atan2(dimag(reflect(nf,ne_previous)),real(reflect(nf,ne_previous), 8))) - - forward_flow(nt)=forward_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& - cos(omega*time+b(nf)+& - atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - - reflected_flow(nt)=reflected_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& - abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& - cos(omega*time+b(nf)+& - atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& - (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& - atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))+& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - - enddo - time=time+dt - enddo - elseif (bc_type.eq.'flow') then - do nt=1,n_time - do nf=1,no_freq - omega=2*pi*nf*harmonic_scale - - forward_pressure(nt)=forward_pressure(nt)+(abs(q_factor(nf,ne))/abs(char_admit(nf,ne)))*a(nf)*& - cos(omega*time+b(nf)+atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))-& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - - forward_pressure_previous(nt)=forward_pressure_previous(nt)+(abs(q_factor(nf,ne_previous))/& - abs(char_admit(nf,ne_previous)))*a(nf)*cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne_previous)),real(q_factor(nf,ne_previous), 8))-& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) + do nt=1,n_time + do nf=1,no_freq + omega=2*pi*nf*harmonic_scale + forward_pressure(nt)=forward_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))) + forward_pressure_previous(nt)=forward_pressure_previous(nt)+abs(p_factor(nf,ne_previous))*& + a(nf)*cos(omega*time+b(nf)+atan2(dimag(p_factor(nf,ne_previous)),real(p_factor(nf,ne_previous), 8))) + + reflected_pressure(nt)=reflected_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*& + abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& + cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& + (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& + atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))) + reflected_pressure_previous(nt)=reflected_pressure_previous(nt)+abs(p_factor(nf,ne_previous))*& + a(nf)*abs(reflect(nf,ne_previous))*exp((-2*elem_field(ne_length,ne_previous))*& + (real(prop_const(nf,ne_previous), 8)))*cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne_previous)),real(p_factor(nf,ne_previous), 8))+& + (-2*elem_field(ne_length,ne_previous))*(dimag(prop_const(nf,ne_previous)))+& + atan2(dimag(reflect(nf,ne_previous)),real(reflect(nf,ne_previous), 8))) + + forward_flow(nt)=forward_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& + cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& + atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) + + reflected_flow(nt)=reflected_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& + abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& + cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& + (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& + atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))+& + atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) + + enddo + time=time+dt + enddo - reflected_pressure(nt)=reflected_pressure(nt)+(abs(q_factor(nf,ne))/abs(char_admit(nf,ne)))*a(nf)*& - abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& - cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))+& - (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& - atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))-& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - - reflected_pressure_previous(nt)=reflected_pressure_previous(nt)+(abs(q_factor(nf,ne_previous))/& - abs(char_admit(nf,ne_previous)))*a(nf)*abs(reflect(nf,ne_previous))*& - exp((-2*elem_field(ne_length,ne_previous))*& - (real(prop_const(nf,ne_previous), 8)))*cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne_previous)),real(q_factor(nf,ne_previous), 8))+& - (-2*elem_field(ne_length,ne_previous))*(dimag(prop_const(nf,ne_previous)))+& - atan2(dimag(reflect(nf,ne_previous)),real(reflect(nf,ne_previous), 8))-& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - - forward_flow(nt)=forward_flow(nt)+abs(q_factor(nf,ne))*a(nf)*& - cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))) - - reflected_flow(nt)=reflected_flow(nt)+abs(q_factor(nf,ne))*a(nf)*& - abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& - cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))+& - (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& - atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))) - - enddo - time=time+dt - enddo - else ! Wrong Boundary condition type - print *, "ERROR: Boundary condition type not recognised. Choose flow or pressure type." - call exit(0) - endif np=elem_nodes(2,ne) ! terminals np_previous=elem_nodes(1,ne) ! upstream node to terminal nodes if(.not.allocated(p_terminal)) allocate (p_terminal(n_time)) @@ -2042,13 +2003,12 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, p_terminal = forward_pressure+reflected_pressure + node_field(nj_bv_press,np) p_previous = forward_pressure_previous+reflected_pressure_previous +& node_field(nj_bv_press,np_previous) - terminal_flow = abs(forward_flow-reflected_flow+elem_field(ne_Qdot,ne)) + terminal_flow = forward_flow - reflected_flow + elem_field(ne_Qdot,ne) + call calculate_ppl(np,grav_vect,mechanics_parameters,Ppl) do nt=1,n_time - delta_p = abs(p_terminal(nt) - p_previous(nt)) - terminals_radius(nt) = sqrt(sqrt((8.0_dp*viscosity*elem_field(ne_length,ne)*terminal_flow(nt))/& - (pi*delta_p))) + delta_p = ((p_terminal(nt) + p_previous(nt))/2) - Ppl + terminals_radius(nt) = elem_field(ne_radius_out0,ne)*(1+alpha(ne)*delta_p) WSS(nt) = (4.0_dp * viscosity * terminal_flow(nt))/(pi * terminals_radius(nt)**3) !wall shear stress at terminals - enddo write(fid,fmt=*) ne, forward_pressure+node_field(nj_bv_press,np) ! incident pressure write(fid2,fmt=*) ne, forward_flow+elem_field(ne_Qdot,ne) ! incident flow @@ -2076,12 +2036,86 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, ne = 1 ! MPA inlet element forward_pressure=0.0_dp ! reseting the for MPA reflected_pressure=0.0_dp ! reseting the for MPA - if (bc_type.eq.'pressure')then + forward_flow=0.0_dp ! reseting the for MPA + reflected_flow=0.0_dp ! reseting the for MPA + elem_rad=0.0_dp + WSS=0.0_dp + do nt=1,n_time + do nf=1,no_freq + omega=2*pi*nf*harmonic_scale + forward_pressure(nt)=forward_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))) + + reflected_pressure(nt)=reflected_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*& + abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& + cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& + (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& + atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))) + + forward_flow(nt)=forward_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& + cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& + atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) + + reflected_flow(nt)=reflected_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& + abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& + cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& + (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& + atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))+& + atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) + + enddo + time=time+dt + enddo + + np=elem_nodes(2,ne) ! end node of an elem + np_previous=elem_nodes(1,ne) ! upstream node of the element + if(.not.allocated(p_terminal)) allocate (p_terminal(n_time)) + if(.not.allocated(p_previous)) allocate (p_previous(n_time)) + if(.not.allocated(elem_flow)) allocate (elem_flow(n_time)) + if(.not.allocated(elem_rad)) allocate (elem_rad(n_time)) + p_terminal = forward_pressure + reflected_pressure + node_field(nj_bv_press,np) + p_previous = node_field(nj_bv_press,np_previous) + elem_flow = forward_flow - reflected_flow + elem_field(ne_Qdot,ne) + call calculate_ppl(np,grav_vect,mechanics_parameters,Ppl) + + do nt=1,n_time + delta_p = ((p_terminal(nt) + p_previous(nt))/2)-Ppl + elem_rad(nt) = elem_field(ne_radius_out0,ne)*(1+alpha(ne)*delta_p) + WSS(nt) = (4.0_dp * viscosity * elem_flow(nt))/(pi * elem_rad(nt)**3) !wall shear stress at terminals + enddo + + write(fid8,fmt=*) forward_pressure+reflected_pressure + node_field(nj_bv_press,np_previous) !Inlet total pressure + write(fid9,fmt=*) forward_flow - reflected_flow + elem_field(ne_Qdot,ne) !Inlet MPA flow + write(fid10,fmt=*) forward_flow !Inlet forward flow + write(fid11,fmt=*) reflected_flow !Inlet reflected flow + write(fid12,fmt=*) forward_pressure !Inlet forward pressure + write(fid13,fmt=*) reflected_pressure !Inlet reflected pressure + write(200,fmt=*) ne, WSS ! WSS for element 1 to be added to WSS_arterial.txt + write(210,fmt=*) ne, elem_rad ! radius for element 1 to be added to radii.txt + write(220,fmt=*) ne, forward_flow - reflected_flow + elem_field(ne_Qdot,ne) ! Flow for element 1 to be added to flows.txt + + + do ne =2,min_ven-1 ! reporting time-dependent radius and WSS for all the arteries + ne_previous=elem_cnct(-1,1,ne) + forward_pressure=0.0_dp + reflected_pressure=0.0_dp + forward_pressure_previous=0.0_dp + reflected_pressure_previous=0.0_dp + forward_flow=0.0_dp + reflected_flow=0.0_dp + WSS=0.0_dp + elem_rad=0.0_dp + elem_flow=0.0_dp do nt=1,n_time do nf=1,no_freq omega=2*pi*nf*harmonic_scale forward_pressure(nt)=forward_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*cos(omega*time+b(nf)+& atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))) + forward_pressure_previous(nt)=forward_pressure_previous(nt)+abs(p_factor(nf,ne_previous))*& + a(nf)*cos(omega*time+b(nf)+atan2(dimag(p_factor(nf,ne_previous)),real(p_factor(nf,ne_previous), 8))) reflected_pressure(nt)=reflected_pressure(nt)+abs(p_factor(nf,ne))*a(nf)*& abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& @@ -2089,6 +2123,12 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, atan2(dimag(p_factor(nf,ne)),real(p_factor(nf,ne), 8))+& (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))) + reflected_pressure_previous(nt)=reflected_pressure_previous(nt)+abs(p_factor(nf,ne_previous))*& + a(nf)*abs(reflect(nf,ne_previous))*exp((-2*elem_field(ne_length,ne_previous))*& + (real(prop_const(nf,ne_previous), 8)))*cos(omega*time+b(nf)+& + atan2(dimag(p_factor(nf,ne_previous)),real(p_factor(nf,ne_previous), 8))+& + (-2*elem_field(ne_length,ne_previous))*(dimag(prop_const(nf,ne_previous)))+& + atan2(dimag(reflect(nf,ne_previous)),real(reflect(nf,ne_previous), 8))) forward_flow(nt)=forward_flow(nt)+abs(char_admit(nf,ne))*abs(p_factor(nf,ne))*a(nf)*& cos(omega*time+b(nf)+& @@ -2106,48 +2146,27 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, enddo time=time+dt enddo - elseif(bc_type.eq.'flow')then - do nt=1,n_time - do nf=1,no_freq - omega=2*pi*nf*harmonic_scale - - forward_pressure(nt)=forward_pressure(nt)+(abs(q_factor(nf,ne))/abs(char_admit(nf,ne)))*a(nf)*& - cos(omega*time+b(nf)+atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))-& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - reflected_pressure(nt)=reflected_pressure(nt)+(abs(q_factor(nf,ne))/abs(char_admit(nf,ne)))*a(nf)*& - abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& - cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))+& - (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& - atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))-& - atan2(dimag(char_admit(nf,ne)),real(char_admit(nf,ne), 8))) - - forward_flow(nt)=forward_flow(nt)+abs(q_factor(nf,ne))*a(nf)*& - cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))) - - reflected_flow(nt)=reflected_flow(nt)+abs(q_factor(nf,ne))*a(nf)*& - abs(reflect(nf,ne))*exp((-2*elem_field(ne_length,ne))*(real(prop_const(nf,ne), 8)))*& - cos(omega*time+b(nf)+& - atan2(dimag(q_factor(nf,ne)),real(q_factor(nf,ne), 8))+& - (-2*elem_field(ne_length,ne))*(dimag(prop_const(nf,ne)))+& - atan2(dimag(reflect(nf,ne)),real(reflect(nf,ne), 8))) - - enddo - time=time+dt + np=elem_nodes(2,ne) ! end node of an elem + np_previous=elem_nodes(1,ne) ! upstream node of the element + if(.not.allocated(p_terminal)) allocate (p_terminal(n_time)) + if(.not.allocated(p_previous)) allocate (p_previous(n_time)) + if(.not.allocated(elem_flow)) allocate (elem_flow(n_time)) + if(.not.allocated(elem_rad)) allocate (elem_rad(n_time)) + p_terminal = forward_pressure + reflected_pressure + node_field(nj_bv_press,np) + p_previous = forward_pressure_previous + reflected_pressure_previous +& + node_field(nj_bv_press,np_previous) + elem_flow = forward_flow - reflected_flow + elem_field(ne_Qdot,ne) ! forwards wave - reflected + S.S flow = total elem flow + call calculate_ppl(np,grav_vect,mechanics_parameters,Ppl) + do nt=1,n_time + delta_p = ((p_terminal(nt) + p_previous(nt))/2) - Ppl + elem_rad(nt) = elem_field(ne_radius_out0,ne)*(1+alpha(ne)*delta_p) + WSS(nt) = (4.0_dp * viscosity * elem_flow(nt))/(pi * elem_rad(nt)**3) !wall shear stress at terminals enddo - else - print *, "ERROR: Boundary condition type not recognised. Choose flow or pressure type." - call exit(0) - endif - np=elem_nodes(1,ne) ! inlet node - write(fid8,fmt=*) forward_pressure+reflected_pressure + node_field(nj_bv_press,np) !Inlet total pressure - write(fid9,fmt=*) forward_flow - reflected_flow + elem_field(ne_Qdot,ne) !Inlet MPA flow - write(fid10,fmt=*) forward_flow !Inlet forward flow - write(fid11,fmt=*) reflected_flow !Inlet reflected flow - write(fid12,fmt=*) forward_pressure !Inlet forward pressure - write(fid13,fmt=*) reflected_pressure !Inlet reflected pressure + write(200,fmt=*) ne, WSS ! elements wall shear stress in cardiac cycle + write(210,fmt=*) ne, elem_rad ! elements radii in cardiac cycle + write(220,fmt=*) ne, forward_flow-reflected_flow+elem_field(ne_Qdot,ne) ! elem flows in cardiac cycle + enddo do ne = 1, num_elems if((elem_field(ne_group,ne).eq.2.0_dp).and.(vein_found.eqv..False.)) then !only applying on veins @@ -2244,6 +2263,14 @@ subroutine evaluate_wave_transmission(grav_dirn,grav_factor,n_time,heartrate,a0, close(fid7) close(fid8) close(fid9) + close(140) + close(150) + close(160) + close(170) + close(180) + close(190) + close(200) + close(210) !!DEALLOCATE MEMORY @@ -2388,12 +2415,13 @@ end subroutine boundary_admittance !############################################################################## ! !*characteristic_admittance* calculates the characteristic admittance of each -subroutine characteristic_admittance(no_freq,char_admit,prop_const,harmonic_scale,& +subroutine characteristic_admittance(no_freq,char_admit,prop_const,alpha,harmonic_scale,& density,viscosity,admit_param,elast_param,mechanics_parameters,grav_vect,remodeling_grade) integer, intent(in) :: no_freq complex(dp), intent(inout) :: char_admit(1:no_freq,num_elems) complex(dp), intent(inout) :: prop_const(1:no_freq,num_elems) + real(dp), intent(inout) :: alpha(num_elems) real(dp), intent(in) :: harmonic_scale real(dp), intent(in) :: density real(dp), intent(in) :: viscosity @@ -2413,6 +2441,7 @@ subroutine characteristic_admittance(no_freq,char_admit,prop_const,harmonic_scal character(len=60) :: sub_name sub_name = 'characteristic_admittance' call enter_exit(sub_name,1) + alpha=elast_param%elasticity_parameters(1) if(remodeling_grade.eq.0.0_dp) then ! Solving for Healthy do ne=1,num_elems do nn=1,2 @@ -2537,8 +2566,10 @@ subroutine characteristic_admittance(no_freq,char_admit,prop_const,harmonic_scal Rg_in=narrow_factor*R0*(Ptm*elast_param%elasticity_parameters(1)+1.d0) elseif(R0.gt.narrow_rad_two) then ! Hypertrophy only Rg_in=R0*(Ptm*alt_hyp*elast_param%elasticity_parameters(1)+1.d0) + alpha(ne)=alt_hyp*elast_param%elasticity_parameters(1) else ! both hypertophy and narrowing Rg_in=narrow_factor*R0*(Ptm*alt_hyp*alt_fib*elast_param%elasticity_parameters(1)+1.d0) + alpha(ne)=alt_hyp*alt_fib*elast_param%elasticity_parameters(1) endif else ! out of range of target vessels,hence, No remodeling Rg_in=R0*(Ptm*elast_param%elasticity_parameters(1)+1.d0) @@ -2550,8 +2581,10 @@ subroutine characteristic_admittance(no_freq,char_admit,prop_const,harmonic_scal Rg_out=narrow_factor*R0*(Ptm*elast_param%elasticity_parameters(1)+1.d0) elseif(R0.gt.narrow_rad_two) then ! only hypertophy Rg_out=R0*(Ptm*alt_hyp*elast_param%elasticity_parameters(1)+1.d0) + alpha(ne)=alt_hyp*elast_param%elasticity_parameters(1) else ! both hypertophy and narrowing Rg_out=narrow_factor*R0*(Ptm*alt_hyp*alt_fib*elast_param%elasticity_parameters(1)+1.d0) + alpha(ne)=alt_hyp*alt_fib*elast_param%elasticity_parameters(1) endif else ! out of range of target vessels,hence, No remodeling Rg_out=R0*(Ptm*elast_param%elasticity_parameters(1)+1.d0) @@ -2819,7 +2852,6 @@ subroutine pressure_flow_factor(no_freq,p_factor,q_factor,reflect,prop_const,cha p_factor=1.0_dp q_factor=1.0_dp - if (bc_type.eq.'pressure') then do nf=1,no_freq omega=nf*2*PI*harmonic_scale do ne=ne_min,ne_max @@ -2827,44 +2859,15 @@ subroutine pressure_flow_factor(no_freq,p_factor,q_factor,reflect,prop_const,cha if(elem_cnct(-1,0,ne).eq.0)then !no upstream elements, inlet, ignore ne_up=ne_min p_factor(nf,ne)=(1.0_dp)!* &!assumes input admittance is the same as characteristic admittance for this vessel - !exp(-1.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))!/& - !(1+reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))) - q_factor(nf,ne)=p_factor(nf,ne)*eff_admit(nf,ne) else ne_up=elem_cnct(-1,1,ne) p_factor(nf,ne)=p_factor(nf,ne_up)*(1+reflect(nf,ne_up))* & exp(-1.0_dp*elem_field(ne_length,ne_up)*prop_const(nf,ne_up))/& (1+reflect(nf,ne)*exp(-2.0_dp*elem_field(ne_length,ne)*prop_const(nf,ne))) - q_factor(nf,ne)=p_factor(nf,ne)*eff_admit(nf,ne) - endif!neup - enddo!ne - enddo!nf - elseif (bc_type.eq.'flow') then ! NEEDS TO BE PROPERLY LOOKED AT. THE MATH ESPECIFICALLY - do nf=1,no_freq - omega=nf*2*PI*harmonic_scale - do ne=ne_min,ne_max - !look for upstream element - if(elem_cnct(-1,0,ne).eq.0)then !no upstream elements, inlet, ignore - ne_up=ne_min - q_factor(nf,ne)=(1.0_dp)!* &!assumes input admittance is the same as characteristic admittance for this vessel - !p_factor(nf,ne)=q_factor(nf,ne)/char_admit(nf,ne)!(1.0_dp)!* &!assumes input admittance is the same as characteristic admittance for this vessel - !exp(-1.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))!/& - !(1+reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))) - else - ne_up=elem_cnct(-1,1,ne) - p_factor(nf,ne)=p_factor(nf,ne_up)*(char_admit(nf,ne_up)*(1-reflect(nf,ne_up))* & - exp(-1.0_dp*elem_field(ne_length,ne_up)*prop_const(nf,ne_up)))/& - (2*char_admit(nf,ne)*(1-reflect(nf,ne)*exp(-2.0_dp*elem_field(ne_length,ne)*prop_const(nf,ne)))) - q_factor(nf,ne)=p_factor(nf,ne_up)*char_admit(nf,ne)!(1+reflect(nf,ne_up))* & - !exp(-1.0_dp*elem_field(ne_length,ne_up)*prop_const(nf,ne_up))/& - !(1+reflect(nf,ne)*exp(-2.0_dp*elem_field(ne_length,ne)*prop_const(nf,ne))) + q_factor(nf,ne)=p_factor(nf,ne)*char_admit(nf,ne) endif!neup enddo!ne enddo!nf - else - print *, "ERROR: Boundary condition not recognised - pressure_flow_factor subroutine" - call exit(0) - endif call enter_exit(sub_name,2) end subroutine pressure_flow_factor From d74b15648491fcc9bcb4d790bea0fc414042c17c Mon Sep 17 00:00:00 2001 From: Behdad Shaarbaf Ebrahimi Date: Tue, 6 Jan 2026 11:52:03 +1300 Subject: [PATCH 2/3] fix bessel function for very small values to not produce NANs --- src/lib/math_utilities.f90 | 272 +++++++++++++++++++++------------- src/lib/wave_transmission.f90 | 4 +- 2 files changed, 167 insertions(+), 109 deletions(-) diff --git a/src/lib/math_utilities.f90 b/src/lib/math_utilities.f90 index 9e9f7a6a..a975919e 100644 --- a/src/lib/math_utilities.f90 +++ b/src/lib/math_utilities.f90 @@ -23,120 +23,178 @@ module math_utilities contains -subroutine bessel_complex(z,bessel0,bessel1) - - complex(dp), intent(in) :: z - complex(dp), intent(out) :: bessel0,bessel1 - - real(dp) :: a(12),a1(10),b(12) - real(dp) :: a0 - complex(dp) :: ci,cr,z1,ca,zr - integer :: k,k0 -! complex ( kind = 8 ) ca -! complex ( kind = 8 ) cb -! complex ( kind = 8 ) ci -! complex ( kind = 8 ) cr -! complex ( kind = 8 ) cs -! complex ( kind = 8 ) ct -! complex ( kind = 8 ) cw -! integer ( kind = 4 ) k -! integer ( kind = 4 ) k0 -! real ( kind = 8 ) pi -! real ( kind = 8 ) w0 -! complex ( kind = 8 ) z -! complex ( kind = 8 ) z1 -! complex ( kind = 8 ) z2 -! complex ( kind = 8 ) zr -! complex ( kind = 8 ) zr2 - a = (/ & - 0.125e00_dp, 7.03125e-02_dp,& - 7.32421875e-02_dp, 1.1215209960938e-01_dp,& - 2.2710800170898e-01_dp, 5.7250142097473e-01_dp,& - 1.7277275025845e00_dp, 6.0740420012735e00_dp,& - 2.4380529699556e01_dp, 1.1001714026925e02_dp,& - 5.5133589612202e02_dp, 3.0380905109224e03_dp /) - a1 = (/ & - 0.125e00_dp, 0.2109375e00_dp, & - 1.0986328125e00_dp, 1.1775970458984e01_dp, & - 2.1461706161499e002_dp, 5.9511522710323e03_dp, & - 2.3347645606175e05_dp, 1.2312234987631e07_dp, & - 8.401390346421e08_dp, 7.2031420482627e10_dp /) - b = (/ & - -0.375e00_dp, -1.171875e-01_dp, & - -1.025390625e-01_dp, -1.4419555664063e-01_dp, & - -2.7757644653320e-01_dp, -6.7659258842468e-01_dp, & - -1.9935317337513e00_dp, -6.8839142681099e00_dp, & - -2.7248827311269e01_dp, -1.2159789187654e02_dp, & - -6.0384407670507e02_dp, -3.3022722944809e03_dp /) + subroutine bessel_complex(z, bessel0, bessel1) + use, intrinsic :: ieee_arithmetic + implicit none -! - ci = cmplx (0.0_dp,1.0_dp,8) - a0 = abs (z) - z1 = z -! - if(abs(a0).le.zero_tol)then - bessel0 = cmplx(1.0_dp,0.0_dp,8) - bessel1 = cmplx(0.0_dp,0.0_dp,8) - endif - - if(real(z).lt.0.0_dp) then - z1 = -z - endif -! - if( a0 <= 18.0_dp ) then - - bessel0 =cmplx(1.0_dp,0.0_dp,8) - cr = cmplx(1.0_dp,0.0_dp,8) - do k = 1,50 - cr = 0.25_dp*cr* z1**2/k**2 - bessel0 = bessel0+cr - if (abs (cr/bessel0).LT.1.0e-15_dp) then - exit - endif - enddo - - bessel1 =cmplx(1.0_dp,0.0_dp,8) - cr = cmplx(1.0_dp,0.0_dp,8) - do k = 1,50 - cr = 0.25_dp*cr*z**2/(k*(k+1)) - bessel1 = bessel1+cr - if (abs (cr/bessel1).LT.1.0e-15_dp) then - exit - endif - enddo - - bessel1 = 0.5_dp*z1*bessel1 - - else - - if ( a0 < 35.0_dp ) then - k0 = 12 - else if ( a0 < 50.0_dp ) then - k0 = 9 + complex(dp), intent(in) :: z + complex(dp), intent(out) :: bessel0, bessel1 + + real(dp) :: a(12), a1(10), b(12) + real(dp) :: a0, absz, scale + complex(dp) :: cr, z1, ca, zr, zwork + integer :: k, k0 + + ! ----------------------- + ! User-tunable safety caps + ! ----------------------- + real(dp), parameter :: Z_SMALL = 1.0e-12_dp ! small-|z| threshold + real(dp), parameter :: Z_ABS_MAX = 200.0_dp ! cap |z| (keeps series/asymptotic stable) + real(dp), parameter :: Z_RE_MAX = 50.0_dp ! cap Re(z) to avoid exp overflow / solver blow-up + real(dp), parameter :: REL_TOL = 1.0e-15_dp + + ! Coefficients (unchanged from your original) + a = (/ & + 0.125e00_dp, 7.03125e-02_dp, & + 7.32421875e-02_dp, 1.1215209960938e-01_dp, & + 2.2710800170898e-01_dp, 5.7250142097473e-01_dp, & + 1.7277275025845e00_dp, 6.0740420012735e00_dp, & + 2.4380529699556e01_dp, 1.1001714026925e02_dp, & + 5.5133589612202e02_dp, 3.0380905109224e03_dp /) + + a1 = (/ & + 0.125e00_dp, 0.2109375e00_dp, & + 1.0986328125e00_dp, 1.1775970458984e01_dp, & + 2.1461706161499e002_dp, 5.9511522710323e03_dp, & + 2.3347645606175e05_dp, 1.2312234987631e07_dp, & + 8.401390346421e08_dp, 7.2031420482627e10_dp /) + + b = (/ & + -0.375e00_dp, -1.171875e-01_dp, & + -1.025390625e-01_dp, -1.4419555664063e-01_dp, & + -2.7757644653320e-01_dp, -6.7659258842468e-01_dp, & + -1.9935317337513e00_dp, -6.8839142681099e00_dp, & + -2.7248827311269e01_dp, -1.2159789187654e02_dp, & + -6.0384407670507e02_dp, -3.3022722944809e03_dp /) + + ! ----------------------- + ! 0) Input sanitisation + ! ----------------------- + if (.not. ieee_is_finite(real(z)) .or. .not. ieee_is_finite(aimag(z))) then + bessel0 = cmplx(0.0_dp, 0.0_dp, kind=dp) + bessel1 = cmplx(0.0_dp, 0.0_dp, kind=dp) + return + end if + + ! ----------------------- + ! 1) Optional clipping + ! (prevents solver-destroying growth from exp(Re(z))) + ! ----------------------- + zwork = z + + ! Clip real part (prevents exp overflow and insane growth) + if (real(zwork) > Z_RE_MAX) zwork = cmplx(Z_RE_MAX, aimag(zwork), kind=dp) + if (real(zwork) < -Z_RE_MAX) zwork = cmplx(-Z_RE_MAX, aimag(zwork), kind=dp) + + ! Clip magnitude (keeps powers zr**k and series stable) + absz = abs(zwork) + if (absz > Z_ABS_MAX .and. absz > 0.0_dp) then + scale = Z_ABS_MAX / absz + zwork = zwork * scale + absz = Z_ABS_MAX + end if + + a0 = absz + + ! ----------------------- + ! 2) Robust small-|z| handling + ! ----------------------- + if (a0 <= max(zero_tol, Z_SMALL)) then + bessel0 = cmplx(1.0_dp, 0.0_dp, kind=dp) ! I0(0)=1 + bessel1 = 0.5_dp * zwork ! I1(z) ~ z/2 near 0 + return + end if + + ! Preserve original symmetry rule based on ORIGINAL z (not clipped zwork) + if (real(z) < 0.0_dp) then + z1 = -zwork else - k0 = 7 + z1 = zwork end if - ca = exp(z1)/sqrt(2.0_dp*pi*z1) - bessel0 = cmplx(1.0_dp,0.0_dp,8) - zr = 1.0_dp/z1 - do k = 1,k0 - bessel0 = bessel0 + a(k) * zr ** k - enddo - bessel0 = ca * bessel0 - bessel1 = cmplx (1.0_dp,0.0_dp,8) - do k = 1,k0 - bessel1 =bessel1+b(k)*zr**k - end do - bessel1 = ca * bessel1 - endif + ! ----------------------- + ! 3) Main evaluation + ! ----------------------- + if (a0 <= 18.0_dp) then + ! ---- Power series (safe region) ---- + + ! I0(z) series + bessel0 = cmplx(1.0_dp, 0.0_dp, kind=dp) + cr = cmplx(1.0_dp, 0.0_dp, kind=dp) + do k = 1, 50 + cr = 0.25_dp * cr * (z1*z1) / real(k*k, dp) + bessel0 = bessel0 + cr + if (abs(cr) < REL_TOL * max(1.0_dp, abs(bessel0))) exit + end do + + ! I1(z) series (match z1 consistently) + bessel1 = cmplx(1.0_dp, 0.0_dp, kind=dp) + cr = cmplx(1.0_dp, 0.0_dp, kind=dp) + do k = 1, 50 + cr = 0.25_dp * cr * (z1*z1) / real(k*(k+1), dp) + bessel1 = bessel1 + cr + if (abs(cr) < REL_TOL * max(1.0_dp, abs(bessel1))) exit + end do + bessel1 = 0.5_dp * z1 * bessel1 + else + ! ---- Asymptotic branch (large |z|) ---- + + if (a0 < 35.0_dp) then + k0 = 12 + else if (a0 < 50.0_dp) then + k0 = 9 + else + k0 = 7 + end if + + ! Guard against tiny z1 (paranoia) + if (abs(z1) <= max(zero_tol, Z_SMALL)) then + bessel0 = cmplx(1.0_dp, 0.0_dp, kind=dp) + bessel1 = 0.5_dp * z1 + return + end if + + ! If exp(z1) would overflow, DO NOT inject huge(). + ! Return bounded values so your network doesn't explode. + if (real(z1) > log(huge(1.0_dp)) - 2.0_dp) then + bessel0 = cmplx(0.0_dp, 0.0_dp, kind=dp) + bessel1 = cmplx(0.0_dp, 0.0_dp, kind=dp) + return + end if + + zr = 1.0_dp / z1 + ca = exp(z1) / sqrt(2.0_dp*pi*z1) + + bessel0 = cmplx(1.0_dp, 0.0_dp, kind=dp) + do k = 1, k0 + bessel0 = bessel0 + a(k) * (zr ** k) + end do + bessel0 = ca * bessel0 + + bessel1 = cmplx(1.0_dp, 0.0_dp, kind=dp) + do k = 1, k0 + bessel1 = bessel1 + b(k) * (zr ** k) + end do + bessel1 = ca * bessel1 + end if - if ( real (z).lt.0.0_dp)then - bessel1 = - bessel1 - endif + ! Preserve original sign correction rule + if (real(z) < 0.0_dp) then + bessel1 = -bessel1 + end if + + ! ----------------------- + ! 4) Final sanitisation (NEVER return NaN/Inf) + ! ----------------------- + if (.not. ieee_is_finite(real(bessel0)) .or. .not. ieee_is_finite(aimag(bessel0))) then + bessel0 = cmplx(0.0_dp, 0.0_dp, kind=dp) + end if + if (.not. ieee_is_finite(real(bessel1)) .or. .not. ieee_is_finite(aimag(bessel1))) then + bessel1 = cmplx(0.0_dp, 0.0_dp, kind=dp) + end if -end subroutine bessel_complex + end subroutine bessel_complex ! diff --git a/src/lib/wave_transmission.f90 b/src/lib/wave_transmission.f90 index e4b536c4..8ce65204 100644 --- a/src/lib/wave_transmission.f90 +++ b/src/lib/wave_transmission.f90 @@ -2710,7 +2710,7 @@ subroutine tree_admittance(no_freq,eff_admit,char_admit,reflect,prop_const,harmo eff_admit(nf,ne)=char_admit(nf,ne)*(1& -reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne)))/& (1+reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))) - else!a terminal + else !a terminal daughter_admit=eff_admit(nf,ne) !a boundary condition is applied here reflect(nf,ne)=(char_admit(nf,ne)-daughter_admit)/& (char_admit(nf,ne)+daughter_admit) @@ -2752,7 +2752,7 @@ subroutine tree_admittance(no_freq,eff_admit,char_admit,reflect,prop_const,harmo eff_admit(nf,ne)=char_admit(nf,ne)*(1& -reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne)))/& (1+reflect(nf,ne)*exp(-2.0_dp*prop_const(nf,ne)*elem_field(ne_length,ne))) - else!a terminal + else !a terminal daughter_admit=eff_admit(nf,ne) !a boundary condition is applied here reflect(nf,ne)=(char_admit(nf,ne)-daughter_admit)/& (char_admit(nf,ne)+daughter_admit) From 267fd14405ed39070999a93613024186f9e5caaa Mon Sep 17 00:00:00 2001 From: Behdad Shaarbaf Ebrahimi Date: Wed, 7 Jan 2026 12:06:48 +1300 Subject: [PATCH 3/3] MacOS Ventura (mac-13) removed from workflow due to being retired --- .github/workflows/build_test.yml | 6 ------ 1 file changed, 6 deletions(-) diff --git a/.github/workflows/build_test.yml b/.github/workflows/build_test.yml index 2e4eb2d9..c2e1875c 100644 --- a/.github/workflows/build_test.yml +++ b/.github/workflows/build_test.yml @@ -34,12 +34,6 @@ jobs: super_lu: true, os: ubuntu-24.04, } - - { - name: "macOS Ventura", - build_type: "Release", - os: macos-13, - gfortran: gfortran-13, - } - { name: "macOS Sonoma", build_type: "Release",