program pipe_expansion
    implicit none
    integer::mat,fixed,iostat_val,i,ns
    double precision::L,Ti,To,OD_mm,wt_mm,E_in
    double precision::alpha,E,Sy,dT,dL,stress,F_anch
    double precision::OD,ID,A_x,I_x,SA,Sc,Sh,L_loop
    double precision::T_sw,dT_sw,dL_sw,stress_sw
    double precision,parameter::pi=3.14159265d0
    character(len=30)::mname,fname
    read(*,*,iostat=iostat_val) mat; read(*,*,iostat=iostat_val) L
    read(*,*,iostat=iostat_val) Ti; read(*,*,iostat=iostat_val) To
    read(*,*,iostat=iostat_val) OD_mm; read(*,*,iostat=iostat_val) wt_mm
    read(*,*,iostat=iostat_val) E_in; read(*,*,iostat=iostat_val) fixed
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    if(L<=0)L=10.0d0;if(OD_mm<=0)OD_mm=114.3d0;if(wt_mm<=0)wt_mm=6.02d0
    select case(mat)
    case(1);mname='Carbon Steel (A106-B)';alpha=12.0d-6;E=200.0d0;Sy=248.0d0
    case(2);mname='304 Stainless Steel';alpha=17.3d-6;E=193.0d0;Sy=205.0d0
    case(3);mname='316 Stainless Steel';alpha=16.0d-6;E=193.0d0;Sy=205.0d0
    case(4);mname='HDPE (PE100)';alpha=150.0d-6;E=0.8d0;Sy=20.0d0
    case default;mname='Copper (C12200)';alpha=17.0d-6;E=117.0d0;Sy=69.0d0;mat=5
    end select
    if(E_in>0)E=E_in
    if(fixed==1)then;fname='Both Ends Fixed (anchored)';else;fname='One End Free (guided)';fixed=2;end if
    dT=To-Ti; OD=OD_mm/1000.0d0; ID=(OD_mm-2.0d0*wt_mm)/1000.0d0
    A_x=pi/4.0d0*(OD**2-ID**2)
    I_x=pi/64.0d0*(OD**4-ID**4)
    dL=alpha*L*dT*1000.0d0  ! mm
    if(fixed==1)then
      stress=E*1000.0d0*alpha*abs(dT)  ! MPa (E in GPa -> *1000)
      F_anch=stress*1.0d6*A_x  ! N
    else
      stress=0.0d0; F_anch=0.0d0
    end if
    Sc=Sy/3.0d0; Sh=Sy/3.0d0
    SA=1.0d0*(1.25d0*Sc+0.25d0*Sh)
    L_loop=0.0d0
    if(abs(dL)>0.0d0.and.SA>0.0d0) &
      L_loop=sqrt(3.0d0*E*1000.0d0*OD*abs(dL)/1000.0d0/SA)
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   PIPE THERMAL EXPANSION & FLEXIBILITY (ASME B31.3)'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')       '  Material                  = ',trim(mname)
    write(*,'(A,F10.2,A)') '  Pipe Length L             = ',L,' m'
    write(*,'(A,F10.2,A)') '  Install Temperature       = ',Ti,' C'
    write(*,'(A,F10.2,A)') '  Operating Temperature     = ',To,' C'
    write(*,'(A,F10.2,A)') '  Pipe OD                   = ',OD_mm,' mm'
    write(*,'(A,F10.2,A)') '  Wall Thickness            = ',wt_mm,' mm'
    write(*,'(A,A)')       '  End Condition             = ',trim(fname)
    write(*,*)
    write(*,'(A)')'--- MATERIAL PROPERTIES -------------------------------------'
    write(*,'(A,ES12.4,A)')'  Expansion Coeff alpha     = ',alpha,' 1/C'
    write(*,'(A,F10.2,A)') '  Elastic Modulus E         = ',E,' GPa'
    write(*,'(A,F10.1,A)') '  Yield Strength Sy         = ',Sy,' MPa'
    write(*,*)
    write(*,'(A)')'--- RESULTS -------------------------------------------------'
    write(*,'(A,F10.2,A)') '  Temperature Difference dT = ',dT,' C'
    write(*,'(A,F10.4,A)') '  Free Expansion dL         = ',dL,' mm'
    write(*,'(A,F10.4,A)') '  Free Expansion dL         = ',dL/25.4d0,' in'
    if(fixed==1)then
    write(*,'(A,F10.2,A)') '  Thermal Stress sigma      = ',stress,' MPa'
    write(*,'(A,F12.1,A)') '  Anchor Force F            = ',F_anch,' N'
    write(*,'(A,F12.1,A)') '  Anchor Force F            = ',F_anch/1000.0d0,' kN'
    else
    write(*,'(A)')         '  Thermal Stress            = 0 (free end)'
    write(*,'(A)')         '  Anchor Force              = 0 (free end)'
    end if
    write(*,'(A,F10.2,A)') '  Allowable Stress Range SA = ',SA,' MPa'
    write(*,'(A,F10.4,A)') '  Expansion Loop Length     = ',L_loop,' m'
    write(*,'(A,ES12.4,A)')'  Cross-Section Area A      = ',A_x,' m2'
    write(*,'(A,ES12.4,A)')'  Moment of Inertia I       = ',I_x,' m4'
    if(fixed==1.and.stress>SA)then
      write(*,'(A)')'  *** WARNING: Stress exceeds allowable range SA ***'
      write(*,'(A)')'  *** Expansion loop or expansion joint required ***'
    end if
    write(*,*)
    ns=30
    write(*,'(A)')'--- EXPANSION VS TEMPERATURE SWEEP --------------------------'
    write(*,'(A)')'  T_oper[C]   dT[C]     dL[mm]      Stress[MPa]'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,ns
      T_sw=Ti+dble(i)*300.0d0/dble(ns)
      dT_sw=T_sw-Ti
      dL_sw=alpha*L*dT_sw*1000.0d0
      if(fixed==1)then;stress_sw=E*1000.0d0*alpha*abs(dT_sw);else;stress_sw=0.0d0;end if
      write(*,'(F8.1,4X,F8.1,4X,F10.4,4X,F10.2)') T_sw,dT_sw,dL_sw,stress_sw
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  dL = alpha * L * dT'
    write(*,'(A)')'  sigma = E * alpha * dT (both fixed)'
    write(*,'(A)')'  SA = f*(1.25*Sc + 0.25*Sh) per ASME B31.3'
    write(*,'(A)')'  L_loop = sqrt(3*E*OD*dL/SA)'
end program pipe_expansion
