c-----------------------------------------------------------------------
c Find steady-state temperature profile across the pipe wall set up by heat
c conduction.
c This program generates a set of matrix coefficients to be piped into a
c general linear equation solver.
c-----------------------------------------------------------------------
c Instructor: Nam Sun Wang
c-----------------------------------------------------------------------

c Declare variables (Divide into 10 subintervals) ----------------------
      parameter (n=10, nmax=100)
      real r(n+1)
      real A(nmax,nmax), b(nmax)

c Output a program header (Nothing, because the output of this is to be piped.)

c Specify the boundary condition ---------------------------------------
      r_in  = 0.05
      r_out = 0.10
      T_in  = 200.
      T_out =  80.

c Call a routine to find the coefficients to linear equations Ax=b -----
      call temp(r_in, r_out, T_in, T_out, n, r, A, nmax, b)

c Output the dimension -------------------------------------------------
      print *, n+1

c Output the coefficients of matrix A & vector b one row at a time -----
      do i=1, n+1
        print *, (A(i,j), j=1, n+1), b(i)
      enddo

      end



c-----------------------------------------------------------------------
      subroutine temp(r_in, r_out, T_in, T_out, n, r, A, lda, b)
c-----------------------------------------------------------------------
c The following differential equation describes the steady-state temperature
c profile across the pipe wall set up by heat conduction.
c
c   dT/dr + r (d^2T/dr^2) = 0       BC: T (r_in)  = T_in
c                                       T (r_out) = T_out
c
c The finite difference representation is:
c
c   (1/2+r/dr)*T(i+1) - 2*r/dr*T(i) + (-1/2+r/dr)*T(i-1) = 0   i=2,3,...,10
c
c where r = r_in+dr*(i-1)       i=1, 2,..., n+1
c
c An example with n=10
c i=1:  T(1) = T_in        ... specified
c i=2:  (-0.5+0.055/0.005)*T(1) -2*0.055/0.005*T(2) + (0.5+0.055/0.005)*T(3) = 0
c i=3:  (-0.5+0.060/0.005)*T(2) -2*0.060/0.005*T(3) + (0.5+0.060/0.005)*T(4) = 0
c i=4:  (-0.5+0.065/0.005)*T(3) -2*0.065/0.005*T(4) + (0.5+0.065/0.005)*T(5) = 0
c  :                                                                         :
c i=10: (-0.5+0.095/0.005)*T(9) -2*0.095/0.005*T(10)+ (0.5+0.095/0.005)*T(11)= 0
c i=11: T(11) = T_out       ... specified
c
c Variables used ...
c   r_in  = radius of the inner wall in meter (input)
c   r_out = radius of the outer wall in meter (input)
c   T_in  = temperature at the inner wall in Celsius (input)
c   T_out = temperature at the outer wall in Celsius (input)
c   n     = number of divisions or intervals (input)
c   r(i)  = radius at point i (output)
c   A     = matrix A as in Ax=b(output)
c   lda   = leading dimension of matrix A (input)
c   b     = vector b as in Ax=b (input)
c-----------------------------------------------------------------------

c Declare variables ----------------------------------------------------
      real r(1)
      real A(lda,1), b(1)

c Assign constants in the linear equation ------------------------------
      dr = (r_out-r_in)/float(n)
      do i=2, n
        r(i)  = r_in + dr*(i-1)
        A(i,i-1) = -1./2. + r(i)/dr
        A(i,i)   = -2.*r(i)/dr
        A(i,i+1) =  1./2. + r(i)/dr
      enddo

c B.C.: at r=r_in, T=T_in ----------------------------------------------
      A(1,1)     = 1.
      b(1)       = T_in

c B.C.: at r=r_out, T=T_out --------------------------------------------
      A(n+1,n+1) = 1.
      b(n+1)     = T_out

      end