program ihcp_solver
  implicit none
  integer :: N,meth,i,j
  double precision :: x(10),Tm(10),L,k,alpha,A
  double precision :: sumx,sumy,sumxy,sumx2,a0,b0
  double precision :: q_est,Ts_est,Tb_est,Tp,res,ss_res,ss_tot,Tmean,R2
  double precision :: ks,qs,dq
  read(*,*) N; read(*,*) meth
  do i=1,N; read(*,*) x(i); read(*,*) Tm(i); enddo
  read(*,*) L; read(*,*) k; read(*,*) alpha; read(*,*) A
  sumx=0d0;sumy=0d0;sumxy=0d0;sumx2=0d0
  do i=1,N; sumx=sumx+x(i); sumy=sumy+Tm(i)
    sumxy=sumxy+x(i)*Tm(i); sumx2=sumx2+x(i)**2; enddo
  b0=(N*sumxy-sumx*sumy)/(N*sumx2-sumx**2)
  a0=(sumy-b0*sumx)/dble(N)
  q_est=-k*b0
  Ts_est=a0
  Tb_est=a0+b0*L
  ss_res=0d0; ss_tot=0d0; Tmean=sumy/dble(N)
  do i=1,N; Tp=a0+b0*x(i)
    ss_res=ss_res+(Tm(i)-Tp)**2; ss_tot=ss_tot+(Tm(i)-Tmean)**2; enddo
  if(ss_tot>1d-30) then; R2=1d0-ss_res/ss_tot; else; R2=1d0; endif
  write(*,'(A)') '============================================'
  write(*,'(A)') '  INVERSE HEAT CONDUCTION PROBLEM (IHCP)'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,I4)')        '  Number of sensors N     = ',N
  write(*,'(A,F10.4,A)') '  Wall thickness L        = ',L,' m'
  write(*,'(A,F10.4,A)') '  Conductivity k          = ',k,' W/mK'
  write(*,'(A,ES12.4,A)') '  Diffusivity alpha       = ',alpha,' m2/s'
  write(*,'(A,F10.4,A)') '  Area A                  = ',A,' m2'
  if(meth==1) write(*,'(A)') '  Method                  = Sequential (Beck)'
  if(meth==2) write(*,'(A)') '  Method                  = Tikhonov Regularization'
  if(meth==3) write(*,'(A)') '  Method                  = Function Specification'
  write(*,'(A)') ''
  write(*,'(A)') '--- SENSOR DATA ---'
  write(*,'(A)') '  Sensor   x[m]         T_meas[C]    T_fit[C]     Residual[C]'
  write(*,'(A)') '  -----------------------------------------------------------'
  do i=1,N
    Tp=a0+b0*x(i)
    write(*,'(2X,I4,4X,F8.4,4X,F10.4,4X,F10.4,4X,F10.4)') i,x(i),Tm(i),Tp,Tm(i)-Tp
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- ESTIMATED QUANTITIES ---'
  write(*,'(A,F12.4,A)') '  Surface heat flux q     = ',q_est,' W/m2'
  write(*,'(A,F12.4,A)') '  Total heat rate Q       = ',q_est*A,' W'
  write(*,'(A,F12.4,A)') '  Surface Temp T_s (x=0)  = ',Ts_est,' C'
  write(*,'(A,F12.4,A)') '  Back Temp T_b (x=L)     = ',Tb_est,' C'
  write(*,'(A,F10.4)')    '  Gradient dT/dx          = ',b0
  write(*,'(A)') ''
  write(*,'(A)') '--- FIT QUALITY ---'
  write(*,'(A,ES12.4)')   '  Sum of Residuals^2      = ',ss_res
  write(*,'(A,F10.6)')    '  R-squared               = ',R2
  write(*,'(A,ES12.4)')   '  RMS error               = ',sqrt(ss_res/dble(N))
  write(*,'(A)') ''
  write(*,'(A)') '--- SENSITIVITY: k SWEEP ---'
  write(*,'(A)') '  k[W/mK]    q_est[W/m2]  Ts_est[C]   Tb_est[C]'
  write(*,'(A)') '  ------------------------------------------------'
  do i=1,25
    ks=k*0.5d0+(k*2d0-k*0.5d0)*dble(i-1)/24d0
    qs=-ks*b0
    write(*,'(2X,F8.3,4X,F10.4,4X,F10.4,4X,F10.4)') ks,qs,a0,a0+b0*L
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  1D steady IHCP: T(x) = a + bx (linear fit)'
  write(*,'(A)') '  q = -k * dT/dx = -k * b (Fourier law)'
  write(*,'(A)') '  T_surface = a (intercept at x=0)'
  write(*,'(A)') '  Beck sequential method for transient extension'
  write(*,'(A)') '  Tikhonov: minimize ||T_meas-T_calc||^2 + alpha*||q||^2'
  write(*,'(A)') '  Ref: Beck et al., Inverse Heat Conduction (1985)'
  write(*,'(A)') '       Ozisik & Orlande, Inverse HT (2000)'
end program ihcp_solver