!  window_temperature - a program to calulate the temperature profile through 
!  a window based on the heat flux, and the conductivities and thicknesses
!  of the materials in the window.
!
!  Uses the pgplot library for plotting 
!
!   By Jim Crumley  -  1/30/2007
!
! Compile with:
! gfortran -lpgplot window_temperature.f90 -o window_temperature

PROGRAM window_temperature
  IMPLICIT NONE
  ! T is temperature in degrees C
  ! const is all of the constants combined ! const = P/A/k W/M^2/ deg C
  ! x is the position of each temperature reading,
  ! delta_x is distance between temperature  locations in m.
  ! P_A is the heat flux in W/m^2
  ! T_i  is the initial temperature in degrees C
  ! num_layers is how many different layers are in the windows
  ! k = conductivity in W/m/degrees C for each layer
  ! width = thickness for each layer
  REAL(KIND=8), DIMENSION(100000)  :: T, const, x
  REAL(KIND=8), PARAMETER :: delta_x=1.d-6, P_A = 400.d0, T_i = -20.
  INTEGER, PARAMETER :: num_layers =3  ! 0.0177 
  REAL(KIND=8), PARAMETER, DIMENSION(num_layers) :: k =(/0.84, 0.0177, 0.84/), &
                                            width =(/0.005, 0.002, 0.005/)
  INTEGER :: i, j,i_delta   ! index variables
  INTEGER, DIMENSION(num_layers) :: i_max, i_min  
  INTEGER :: ier       ! check for error with plot file + pgplot function
  INTEGER, EXTERNAL :: pgbeg

     ! minimimum and maximum values for index in each region
 
  !find indices binding each region.
  i_min(1) = 1
  DO i=1,num_layers
    i_max(i) = SUM(width(1:i))/ delta_x
  END DO
  i_min(2:num_layers) = i_max(1:num_layers-1) - 1

  
  !setup x for grid
  DO i = 1, i_max(num_layers)
    x(i) = (i) * delta_x 
  END DO

  ! setup const
  DO i = 1, num_layers
    const(i_min(i):i_max(i)) = P_A / k(i)
  END DO


  !set temperature at the outside of the window
  T(1) = T_i

  ! Loop through and find temperature elsewhere
  DO i = 2, i_max(num_layers)
    T(i) = T(i-1) + const(i)*delta_x
  ENDDO
  
  PRINT *, 'Temperatures at the edge of each layer ', T(i_max(1)), &
    T(i_max(2)), T(i_max(num_layers))



  ! Plotting section - note that bgplot needs kind =4 REALs - use the
  ! REAL function to convert higher precision numbers to the correct type

  ! Initialize plot
  ier = pgbeg(0, '/xserve',1,1) !
  !ier= pgbeg(0,'?',1,1)

  ! set scale for plot
  CALL pgenv(0., MAXVAL(REAL(x(1:i_max(num_layers)))), REAL(T_i),  &
    MAXVAL(REAL(T(1:i_max(num_layers)))), 0, 0)

  !set labels
  CALL pglab ('position(m)', 'T (degrees C)', 'Temperature profile for window')

  !plot the points
  CALL pgpt(i_max(num_layers), REAL(x(1:i_max(num_layers))), &
    REAL(T(1:i_max(num_layers))), 1)
  
  !close the plot
  CALL pgend

END PROGRAM
