;HW5: Factorial and Linear Interpolation
function fac,x
   if x eq 0 then return,1 else return,float(x)*fac(x-1)
end

pro linterp,x,y,z,x2,y2,z2
   nx=n_elements(x)
   ny=n_elements(y)
   x2=fltarr(nx*2-1)
   y2=fltarr(ny*2-1)
   z2=fltarr(nx*2-1,ny*2-1)
   for j=0,ny-1 do begin
      y2[j*2]=y[j]
      if 2*j+1 lt ny*2-1 then y2[j*2+1]=(y[j]+y[j+1])/2.
   endfor
   for j=0,nx-1 do begin
      for l=0,ny-1 do begin
	z2[j*2,l*2]=z[j,l]
	if 2*j+1 lt nx*2-1 then z2[j*2+1,l*2]=(z[j,l]+z[j+1,l])/2.
	if 2*l+1 lt ny*2-1 then z2[j*2,l*2+1]=(z[j,l]+z[j,l+1])/2.
	if 2*j+1 lt nx*2-1 AND 2*l+1 lt ny*2-1 then z2[j*2+1,l*2+1]=(z[j,l]+z[j+1,l]+z[j,l+1]+z[j+1,l+1])/4.
      endfor
      x2[j*2]=x[j]
      if 2*j+1 lt nx*2-1 then x2[j*2+1]=(x[j]+x[j+1])/2.
   endfor
end

for j=0,10 do print,fac(j)
x=findgen(5)
y=findgen(5)
z=fltarr(5,5)
openr,1,'contour.dat'
readf,1,z
contour,z,x,y,/follow
linterp,x,y,z,x2,y2,z2
window,1
contour,z2,x2,y2,/follow
linterp,x2,y2,z2,x4,y4,z4
window,2
contour,z4,x4,y4,/follow
end