MODULE Power_Cost_Functions
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Network_Data_Structures
      USE Power_Cost_Parameters
      IMPLICIT NONE
      PUBLIC::PowerElementalCosts
      PRIVATE
      
      CONTAINS
      
      
      SUBROUTINE PowerElementalCosts(flow,potential,resistance,cost,&
      arc_index,arguments_status,tolerance)
      IMPLICIT NONE
      REAL(KIND=r_wp),INTENT(INOUT)::flow,potential,resistance,cost
      INTEGER(KIND=i_wp),INTENT(IN)::arc_index
      CHARACTER(LEN=4),INTENT(IN)::arguments_status
      REAL(KIND=r_wp),INTENT(IN)::tolerance
      
      LOGICAL::known_flow,known_potential
      CHARACTER::status
      REAL(KIND=r_wp)::eps
      
      eps=MAX(10.0*EPSILON(1.0_r_wp),tolerance)
      
      
      status=arguments_status(1:1)
      known_flow=(status=='F')
      IF((.NOT.known_flow).AND.status/='D')THEN
      flow=SIGN(1.0_r_wp,potential)*&
      (ABS(potential)/(alpha*arcs_cost_parameters(1,arc_index)+eps))**(inv_alpha_minus_1)
      
      known_flow=.TRUE.
      END IF
      
      
      status=arguments_status(2:2)
      known_potential=(status=='F')
      IF((.NOT.known_potential).AND.status/='D')THEN
      potential=alpha*arcs_cost_parameters(1,arc_index)*&
      SIGN(1.0_r_wp,flow)*((ABS(flow)+eps)**(alpha-1.0_r_wp))
      
      known_potential=.TRUE.
      END IF
      
      
      status=arguments_status(3:3)
      IF(status/='F'.OR.status/='D')THEN
      IF(known_potential.AND.known_flow)THEN
      resistance=(alpha-1.0_r_wp)*ABS(potential)/(ABS(flow)+eps)
      ELSE IF(known_potential)THEN
      resistance=(alpha-1.0_r_wp)*ABS(potential)*&
      (ABS(potential)/(alpha*arcs_cost_parameters(1,arc_index)+eps)+eps)**(-inv_alpha_minus_1)
      ELSE IF(known_flow)THEN
      resistance=alpha*(alpha-1.0_r_wp)*&
      arcs_cost_parameters(1,arc_index)*((ABS(flow)+eps)**(alpha-2.0_r_wp))
      ELSE
      resistance=0.0_r_wp
      END IF
      END IF
      
      
      status=arguments_status(4:4)
      IF(status/='F'.OR.status/='D')THEN
      IF(known_potential.AND.known_flow)THEN
      cost=(potential*flow)/alpha
      ELSE IF(known_potential)THEN
      cost=arcs_cost_parameters(1,arc_index)*&
      (ABS(potential)/(alpha*arcs_cost_parameters(1,arc_index)+eps)+eps)**&
      (alpha*inv_alpha_minus_1)
      ELSE IF(known_flow)THEN
      cost=arcs_cost_parameters(1,arc_index)*ABS(flow)**alpha
      ELSE
      cost=0.0_r_wp
      END IF
      END IF
      
      END SUBROUTINE PowerElementalCosts
      
      
      
      
      
      
      
      
      
      
      
      END MODULE Power_Cost_Functions