module charge_rpn_idon !rpn reverse polish notation
implicit none
integer,parameter,private :: DB_out_dim=17
character(len=7),dimension(DB_out_dim),private :: DB_out=(/&
'posZ ','time ','En ','sigZ ','DEn ','emitZ ','zXEpMed',&
'Xmed ','sigX ','divergX','emitX ','xXxpMed',&
'Ymed ','sigY ','divergY','emitY ','yXypMed'/)
real,pointer,dimension(:),private :: P_out_parz,P_out_all
!character(len=39) :: var_string='sigx sigy sigz emitx emity emitz En Den'
type Z_Dynamics
real :: pos_z,time,Energy,sig,DEnergy,emit,zXEp_med
end type
type X_Y_Dynamics
real :: pos_z,time,med ,sig,diverg ,emit,xXxp_med
end type
real,allocatable,private :: stack(:),stackcpy(:)
integer,private :: stacksize
!real, pointer,dimension(1:21) ::P_variable !3x7=21 Max. Numb. of var. in Astra Output
contains
subroutine load_idon_rnp_equation(input_file,idon_equation)
implicit none
character(len=12) :: input_file
character(len=10) :: key_name='[idoneity]'
character(len=80):: idon_equation,line
integer :: i,j,file_record,stat_200,N_record
N_record=file_record(input_file)
open (200,file=input_file,iostat=stat_200)
j=0
do i=1,N_record
read (200,fmt='(A)',iostat=stat_200) line
if (stat_200/=0) exit
if (index(line,key_name)==0) cycle
idon_equation=line
read (200,'(A255)') idon_equation
do while (trim(adjustl(idon_equation))=="")
read (200,'(A255)') idon_equation
enddo
print *,'rpn idno equation=',idon_equation
exit
end do
close (200)
if (index(line,key_name)==0) then
print *,'ATTENTION: No Key name [idoneity] into the GIN.in file'
stop
end if
end subroutine
function idon_evaluation(Xout,Yout,Zout,idon_equation) result(idon_val)
implicit none
type (X_Y_Dynamics),target :: Xout,Yout
type (Z_Dynamics) ,target :: Zout
character(len=80),intent(in) :: idon_equation
character(len=80) :: input
integer :: cnt,i,j
real :: r,r1,r2,idon_val
input=idon_equation
allocate (P_out_parz(1:DB_out_dim),P_out_all(1:DB_out_dim))
P_out_parz(1)=Zout%pos_z; P_out_parz(2)=Zout%time; P_out_parz(3)=Zout%Energy
P_out_parz(4)=Zout%sig; P_out_parz(5)=Zout%Denergy; P_out_parz(6)=Zout%emit
P_out_parz(7)=Zout%zXEp_med
P_out_parz(8)=Xout%med; P_out_parz(9)=Xout%sig; P_out_parz(10)=Xout%diverg
P_out_parz(11)=Xout%emit; P_out_parz(12)=Xout%xXxp_med
P_out_parz(13)=Yout%med; P_out_parz(14)=Yout%sig; P_out_parz(15)=Yout%diverg
P_out_parz(16)=Yout%emit; P_out_parz(17)=Yout%xXxp_med
P_out_all=>P_out_parz
stacksize=0
cnt=0 !! initialize, so we can do a +1 the fist time around
outer: do
! truncate beginning of input line, so we can cycle to the next space.
input=trim(adjustl(input(cnt+1:)))
if (len_trim(input) <= 0) then ! exit if input line ends.
exit
end if
cnt=index(input, ' ') ! search for first space in input line
!---search for variables---
inner: do i=1,DB_out_dim
if (index('0123456789+-/*',input(1:1))>=1) exit inner !se è un numero mi butta subito fuori
if (index('expsqr',trim(input(:cnt)))>=1) exit inner !se è exp o sqr mi butta fuori dal ciclo
!print *,'inner_do ',trim(DB_out(i)),' ',input(:cnt),P_out_all(i)
if (trim(DB_out(i))==input(:cnt)) then
r=P_out_all(i)
call push(r)
cycle outer
end if
end do inner
!---search for operators---
select case(input(:cnt))
case('sqr')
call pop(r)
r=r**2
call push(r)
case('exp')
call pop(r)
r=exp(r)
call push(r)
case('-')
call pop(r1)
call pop(r2)
r = r2 - r1
call push(r)
cycle outer
case('+')
call pop(r1)
call pop(r2)
r = r2 + r1
call push(r)
cycle outer
case('*')
call pop(r1)
call pop(r2)
r = r2 * r1
call push(r)
cycle outer
case('/')
call pop(r1)
call pop(r2)
if (r1 == 0) then
write(*,*)'AIIEEE! Division by zero!'
stop
end if
r = r2 / r1
call push(r)
cycle outer
!---search for numbers---
case default
!! okay, we have an number, push to stack
!print *,'rpn defalut=',input(:cnt)
read(input(:cnt),*)r
call push(r)
cycle outer
end select
!! all valid cases should be considered above !!
!write(*,*)'Invalid character in input'
end do outer
idon_val=stack(stacksize)
!print *,'idon_val=',idon_val
end function
subroutine push(r) ! push r ontop of stack
implicit none
real,intent(in) :: r
if (allocated(stack)) then
allocate (stackcpy(stacksize))
stackcpy(:)=stack(:)
deallocate (stack)
stacksize=stacksize+1
allocate (stack(stacksize))
stack(1:(stacksize-1))=stackcpy(:)
stack(stacksize) =r
deallocate (stackcpy)
else
stacksize=stacksize+1
allocate (stack(stacksize))
stack(stacksize)=r
end if
end subroutine
subroutine pop(r) ! pop r from stack
implicit none
real,intent(out) :: r
r=stack(stacksize)
stacksize=stacksize-1
if (allocated(stackcpy)) deallocate (stackcpy)
allocate (stackcpy(stacksize))
stackcpy(:)=stack(:stacksize)
deallocate (stack)
allocate (stack(stacksize))
stack=stackcpy
deallocate (stackcpy)
end subroutine
end module
