PROGRAM IV_xi
      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
      
      INTEGER,PARAMETER::n_xi_points=7
      CHARACTER(LEN=25)::legend=""
      REAL(KIND=r_wp)::value,xi,interval(2),xi_interval(2)=(/0.001_r_wp,0.25_r_wp/),axis(4)
      REAL(KIND=r_wp),DIMENSION(:),ALLOCATABLE::flow,voltage,resistance
      INTEGER::n_points,point,xi_point
      CHARACTER(LEN=4)::arguments_status
      LOGICAL::change_flow
      
      CALL StartProgram
      CALL InitializeLatticeNetworkProblem
      CALL CreateLatticeNetworkProblem(create_tree=.TRUE.)
      
      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
      
      CALL InitGraphics(file="VI_xi.ps",file_type="CONS",&
      n_plots=n_xi_points,tick_labels="4F2",legend_position="UL",&
      plot_title=(/"Flow-Tension Characterstics with Critical Behaviour",&
      "(for several regularization parameters $\xi$)"/),&
      x_label="flow $x$ (tension $t$)",y_label="tension $t$ (flow $x$)")
      
      XiLoop:DO xi_point=1,n_xi_points
      
      xi=xi_interval(1)+(xi_point-1)*(xi_interval(2)-xi_interval(1))/(n_xi_points-1)
      width_distribution='Fixed'
      width_parameters=(/xi,0.0_r_wp/)
      CALL AssignCostParameters
      
      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
      
      
      IF(xi_point==1)axis=(/interval(1),interval(2),0.0_r_wp,MAXVAL(voltage)/)
      WRITE(legend,"(A,F6.3,A)")"$\xi=",xi,"$"
      CALL Plot2D(x=flow,y=voltage,line_spec=(/xi_point,MAX(INT(2.0*n_xi_points/xi_point),1)/),&
      color_spec=(/0,xi_point/),axis=axis,legend=TRIM(legend))
      
      END DO XiLoop
      CALL EndGraphics()
      
      
      END PROGRAM IV_xi