PROGRAM IV_test
      USE Precision
      USE Error_Handling
      USE System_Monitors
      USE Initialization_Termination
      USE Network_Data_Structures
      USE Lattice_Geometry
      USE Network_Geometry
      USE Simple_Graphics
      USE Network_Graphics
      USE Network_Spanning_Trees
      USE jj_i_Cost_Parameters
      USE jj_i_Cost_Functions
      USE Lattice_Network_Optimization
      USE Conjugate_Gradient
      USE Dual_Network_Solvers
      USE Dual_Line_Minimizers
      USE Dual_Newton_SSCNO
      USE SSCNO_Interface
      
      IMPLICIT NONE
      
      REAL(KIND=r_wp)::value,interval(2)
      REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::flow,voltage,resistance
      INTEGER::n_points,point
      CHARACTER(LEN=4)::arguments_status
      LOGICAL::change_flow
      
      CALL StartProgram
      CALL InitializeLatticeNetworkProblem
      CALL CreateLatticeNetworkProblem(create_tree=.TRUE.)
      CALL AssignCostParameters
      
      WRITE(*,*)"Change flow (or change_voltage)?"
      READ(*,*)change_flow
      WRITE(*,*)"Enter desired interval:"
      READ(*,*)interval
      WRITE(*,*)"Enter the number of points"
      READ(*,*)n_points
      
      ALLOCATE(flow(n_points),voltage(n_points),resistance(n_points))
      
      IF(change_flow)THEN
      arguments_status="FCCD"
      ELSE
      arguments_status="CFCD"
      END IF
      
      DO point=1,n_points
      value=interval(1)+REAL(point-1)*(interval(2)-interval(1))/REAL(n_points-1)
      IF(change_flow)THEN
      flow(point:point)=value
      ELSE
      voltage(point:point)=value
      END IF
      CALL jj_iElementalCosts(&
      arguments_status=arguments_status,arcs_flows=flow(point:point),&
      arcs_voltages=voltage(point:point),arcs_resistances=resistance(point:point),&
      arcs_indices=(/1,1/),tolerance=EPSILON(1.0_r_wp))
      END DO
      
      WRITE(*,"(3A30)")"I","V","R"
      WRITE(*,*)"_________________________________________________________________"
      DO point=1,n_points
      WRITE(*,"(3F30.16)")flow(point),voltage(point),resistance(point)
      END DO
      WRITE(*,*)"_________________________________________________________________"
      
      CALL InitGraphics(file="IV.png",file_type="CONS",&
      n_plots=1,tick_labels="9F1",legend_position="UL",&
      plot_title=(/"V(I)"/),&
      x_label="$I$",y_label="$V$")
      CALL Plot2D(x=flow,y=voltage,plot_spec="SCR")
      CALL Plot2D(x=flow,y=voltage,plot_spec="L-B")
      CALL EndGraphics()
      
      CALL InitGraphics(file="IV.png",file_type="CONS",&
      n_plots=1,tick_labels="9F1",legend_position="UL",&
      plot_title=(/"R(I)"/),&
      x_label="$I$",y_label="$R$")
      CALL Plot2D(x=flow,y=resistance,plot_spec="SCR")
      CALL Plot2D(x=flow,y=resistance,plot_spec="L-B")
      CALL EndGraphics()
      
      CALL InitGraphics(file="IV.png",file_type="CONS",&
      n_plots=1,tick_labels="9F1",legend_position="UL",&
      plot_title=(/"R(V)"/),&
      x_label="$V$",y_label="$R$")
      CALL Plot2D(x=voltage,y=resistance,plot_spec="SCR")
      CALL Plot2D(x=voltage,y=resistance,plot_spec="L-B")
      CALL EndGraphics()
      
      END PROGRAM IV_test