goto,weiter
  ;;; read frames
f=findfile('../FITS/RedCont/SunScan-rc*', count=num)
pp=fltarr(1238, 991, num)
ss=intarr(2, num)
for i=0,num-1 do begin
  rdfits,p,f(i)
  p=p(*,*,2)
  ss(*,i)=size(p,/dim)-1
  pp(0:ss(0,i),0:ss(1,i),i)=p
end

stop

  ;;; get the gradient (estimate)
grad=total(pp(0:1200,0:936,2:36),3)

stop

ssh=intarr(2,num)
for i=1,num-1 do begin
  show,pp(*,*,i),2,sc=.5,/use_ex
  show,pp(*,*,i-1),sc=.5,/use_ex
  print,"Mark structure in frame 1"
  cursor,x1,y1,/dev
  wait, .5
  print,"Mark same structure in frame 2"
  wset,2
  cursor,x2,y2,/dev
  dy=2*(y2-y1)
  dx=2*(x2-x1)
  ssh(*,i)=-shc(pp((-dx>0):ss(0,i-1),0:ss(1,i-1),i-1), $
                pp((dx>0):ss(0,i),dy:ss(1,i),i),/filt) + [0, dy]
  print,i,": ",ssh(*,i)
end

for i=1,num-1 do ssh(*,i)=ssh(*,i-1)-ssh(*,i)
ssh(0,*)-=min(ssh(0,*))
ssh(1,*)-=min(ssh(1,*))

weiter:

big=fltarr(max(ss(0,*)+ssh(0,*)+1), ssh(1,0)+ss(1,0)+1)
w=fltarr(max(ss(0,*)+ssh(0,*)+1), ssh(1,0)+ss(1,0)+1)

sub_fit=1
sc=2.5
  ;;; slopes:
cx=[1.,(2.6129117e-06)/sc]
cy=[1.,(4.5194362e-06)/sc]

im=19
p1=pp(0:ss(0,im),0:ss(1,im),im)
if sub_fit eq 1 then begin
    fit=poly(indgen(ss(0,im)+1)-ss(0,im)/2,cx) $
        # poly(indgen(ss(1,im)+1)-ss(1,im)/2, cy)
    p1/=fit
endif

paste,big,p1,w,ssh(0,im),ssh(1,im)

norm=dblarr(num)
norm(19)=1.
cpar=[1.,(4.1589181e-05)/1.6]
for im=18,0,-1 do begin
      ;;; the 4 corners of the overlap
    x1=max(ssh(0,im:im+1))
    x2=min(ssh(0,im:im+1)+ss(0,im:im+1))
    y1=ssh(1,im)
    y2=ssh(1,im+1)+ss(1,im+1)
      ;;; image coordinates im
    xa1=x1-ssh(0,im)
    xa2=x2-ssh(0,im)
    ya1=0
    ya2=y2-ssh(1,im)
      ;;; image coordinates ref
    xb1=x1-ssh(0,im+1)
    xb2=x2-ssh(0,im+1)
    yb1=y1-ssh(1,im+1)
    yb2=ss(1,im+1)

    p1=pp(0:ss(0,im),0:ss(1,im),im)
    p2=pp(0:ss(0,im+1),0:ss(1,im+1),im+1)
    if sub_fit eq 1 then begin
        fit=poly(indgen(ss(0,im)+1)-ss(0,im)/2, cx) $
            # poly(indgen(ss(1,im)+1)-ss(1,im)/2, cy)
        p1/=fit
        fit=poly(indgen(ss(0,im+1)+1)-ss(0,im+1)/2, cx) $
            # poly(indgen(ss(1,im+1)+1)-ss(1,im+1)/2, cy)
        p2/=fit
    endif

    i1=avg(p1(xa1:xa2,ya1:ya2))
    i2=avg(p2(xb1:xb2,yb1:yb2))
    norm(im)=norm(im+1)*i2/i1
    
    paste,big,p1*norm(im),w,ssh(0,im),ssh(1,im)

end

for im=20,num-1 do begin
      ;;; the 4 corners of the overlap
    x1=max(ssh(0,im-1:im))
    x2=min(ssh(0,im-1:im)+ss(0,im-1:im))
    y1=ssh(1,im-1)
    y2=ssh(1,im)+ss(1,im)
      ;;; image coordinates im
    xa1=x1-ssh(0,im)
    xa2=x2-ssh(0,im)
    ya1=y1-ssh(1,im)
    ya2=y2-ssh(1,im)
      ;;; image coordinates ref
    xb1=x1-ssh(0,im-1)
    xb2=x2-ssh(0,im-1)
    yb1=y1-ssh(1,im-1)
    yb2=y2-ssh(1,im-1)

    p1=pp(0:ss(0,im),0:ss(1,im),im)
    p2=pp(0:ss(0,im-1),0:ss(1,im-1),im-1)
    if sub_fit eq 1 then begin
        fit=poly(indgen(ss(0,im)+1)-ss(0,im)/2, cx) $
            # poly(indgen(ss(1,im)+1)-ss(1,im)/2, cy)
        p1/=fit
        fit=poly(indgen(ss(0,im-1)+1)-ss(0,im-1)/2, cx) $
            # poly(indgen(ss(1,im-1)+1)-ss(1,im-1)/2, cy)
        p2/=fit
    endif

    i1=avg(p1(xa1:xa2,ya1:ya2))
    i2=avg(p2(xb1:xb2,yb1:yb2))
    norm(im)=norm(im-1)*i2/i1
    
    paste,big,p1*norm(im),w,ssh(0,im),ssh(1,im)

end

  ;;; rim in top and bottom image:  
rt=794
rb=325
yy=ssh(1,*)+ss(1,*)/2

plot,yy-rb,norm & oplot,ssh(1,0)+rt-yy,norm,/line

end