program anisotropic_cond
  implicit none
  integer :: i
  double precision :: kx,ky,kz,Lx,Ly,Lz,Th,Tc,Ac,theta
  double precision :: dT,qx,qy,qz,Qx_rate,Qy_rate,Qz_rate,Rx,Ry,Rz
  double precision :: keff,th_rad,ths,keffs,pi
  pi=3.14159265358979d0
  read(*,*) kx; read(*,*) ky; read(*,*) kz
  read(*,*) Lx; read(*,*) Ly; read(*,*) Lz
  read(*,*) Th; read(*,*) Tc; read(*,*) Ac; read(*,*) theta
  dT=abs(Th-Tc)
  th_rad=theta*pi/180d0
  qx=kx*dT/Lx; qy=ky*dT/Ly; qz=kz*dT/Lz
  Qx_rate=qx*Ac; Qy_rate=qy*Ac; Qz_rate=qz*Ac
  Rx=Lx/(kx*Ac); Ry=Ly/(ky*Ac); Rz=Lz/(kz*Ac)
  keff=kx*cos(th_rad)**2+ky*sin(th_rad)**2
  write(*,'(A)') '============================================'
  write(*,'(A)') '  ANISOTROPIC THERMAL CONDUCTION'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,F12.4,A)') '  k_x                     = ',kx,' W/mK'
  write(*,'(A,F12.4,A)') '  k_y                     = ',ky,' W/mK'
  write(*,'(A,F12.4,A)') '  k_z                     = ',kz,' W/mK'
  write(*,'(A,F10.4,A)') '  L_x                     = ',Lx,' m'
  write(*,'(A,F10.4,A)') '  L_y                     = ',Ly,' m'
  write(*,'(A,F10.4,A)') '  L_z                     = ',Lz,' m'
  write(*,'(A,F10.2,A)') '  T_hot                   = ',Th,' C'
  write(*,'(A,F10.2,A)') '  T_cold                  = ',Tc,' C'
  write(*,'(A,F10.4,A)') '  Cross Section A         = ',Ac,' m2'
  write(*,'(A,F10.2,A)') '  Rotation theta          = ',theta,' deg'
  write(*,'(A)') ''
  write(*,'(A)') '--- DIRECTIONAL PROPERTIES ---'
  write(*,'(A,F12.4,A)') '  Heat flux q_x           = ',qx,' W/m2'
  write(*,'(A,F12.4,A)') '  Heat flux q_y           = ',qy,' W/m2'
  write(*,'(A,F12.4,A)') '  Heat flux q_z           = ',qz,' W/m2'
  write(*,'(A,F12.4,A)') '  Heat rate Q_x           = ',Qx_rate,' W'
  write(*,'(A,F12.4,A)') '  Heat rate Q_y           = ',Qy_rate,' W'
  write(*,'(A,F12.4,A)') '  Heat rate Q_z           = ',Qz_rate,' W'
  write(*,'(A,F12.6,A)') '  Resistance R_x          = ',Rx,' K/W'
  write(*,'(A,F12.6,A)') '  Resistance R_y          = ',Ry,' K/W'
  write(*,'(A,F12.6,A)') '  Resistance R_z          = ',Rz,' K/W'
  write(*,'(A)') ''
  write(*,'(A)') '--- ANISOTROPY RATIOS ---'
  write(*,'(A,F10.4)')    '  k_x / k_y               = ',kx/ky
  write(*,'(A,F10.4)')    '  k_x / k_z               = ',kx/kz
  write(*,'(A,F10.4)')    '  k_y / k_z               = ',ky/kz
  write(*,'(A,F12.4,A)') '  k_eff(theta) in x-y     = ',keff,' W/mK'
  write(*,'(A)') ''
  write(*,'(A)') '--- ROTATED k_eff (x-y PLANE) ---'
  write(*,'(A)') '  Tensor rotation: k_eff = kx*cos^2(theta) + ky*sin^2(theta)'
  write(*,'(A)') ''
  write(*,'(A)') '--- ANGLE SWEEP ---'
  write(*,'(A)') '  theta[deg]  k_eff[W/mK]   q[W/m2]     Q[W]'
  write(*,'(A)') '  -----------------------------------------------'
  do i=1,19
    ths=dble(i-1)*5d0
    keffs=kx*cos(ths*pi/180d0)**2+ky*sin(ths*pi/180d0)**2
    write(*,'(2X,F8.1,4X,F10.4,4X,F10.4,4X,F10.4)') ths,keffs,keffs*dT/Lx,keffs*dT/Lx*Ac
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  k_eff(theta) = kx*cos^2(theta) + ky*sin^2(theta)'
  write(*,'(A)') '  Off-diagonal: k_xy = (kx-ky)*sin(theta)*cos(theta)'
  write(*,'(A)') '  3D rotation requires full tensor transformation'
  write(*,'(A)') '  Ref: Incropera Appendix A, Ozisik Ch.1'
  write(*,'(A)') '       Chung, Thermal Conductivity of Composites'
end program anisotropic_cond