subroutine align
c
c	determine les parametres des equations de passage entre les 
c	systeme cartesien de l'afficheur et de gsc
c
	include 	'../common.cod'
	include 	'../common.par'
	include 	'../erreur.par'
	include 	'../iblock.cod'
	include		'align.key'
	logical		flag
	integer		ilen
	character*80	message
	character*20	name

	parameter	(ns_max=20)
	integer*4	stat,mark(ns_max)
	real*4		xres(ns_max),yres(ns_max)
	real*8		tr(6),f1,f2,ang(2),rms
        real*4		rms0
	real*8		pi

        data pi/3.1415927/
	
	o_alif1  = i_alif1
	o_alif2  = i_alif2

	o_alisig = i_alisig
	call ifqnnn('FSIG',1,o_alisig,flag)

	o_alimet = i_alimet
        call ifqcnn('METHODE',1,o_alimet,ilen,flag)
	call maj_min(o_alimet)

	call ifpara(1,message,ilen,flag)
	read(message,*,iostat=ios)ns
	if(ios.ne.0 .or. ns.le.0 .or. ns.gt.ns_max)then
	  write(*,'(1x,a,a)')'Nombre de coordonnees illegale dans: ',
     .               message(1:ilen)
	  goto 9999
	endif

	name = 'X1'
        call ifqcnn('X1',1,name,ilen,flag)
	call read_num_kw(name,1,qq,iad_x1)
	if(erreur)call errare(e_nokblo,6,name,' ',1,0,*9999)
        call test_lim(name,'N',1,ns)
	if(erreur)call errare(e_kwnoli,6,name,' ',1,0,*9999)
	
	name = 'Y1'
        call ifqcnn('Y1',1,name,ilen,flag)
	call read_num_kw(name,1,qq,iad_y1)
	if(erreur)call errare(e_nokblo,6,name,' ',1,0,*9999)
        call test_lim(name,'N',1,ns)
	if(erreur)call errare(e_kwnoli,6,name,' ',1,0,*9999)
	
	name = 'X2'
        call ifqcnn('X2',1,name,ilen,flag)
	call read_num_kw(name,1,qq,iad_x2)
	if(erreur)call errare(e_nokblo,6,name,' ',1,0,*9999)
        call test_lim(name,'N',1,ns)
	if(erreur)call errare(e_kwnoli,6,name,' ',1,0,*9999)
	
	name = 'Y2'
        call ifqcnn('Y2',1,name,ilen,flag)
	call read_num_kw(name,1,qq,iad_y2)
	if(erreur)call errare(e_nokblo,6,name,' ',1,0,*9999)
        call test_lim(name,'N',1,ns)
	if(erreur)call errare(e_kwnoli,6,name,' ',1,0,*9999)

        f1 = i_alif1
	
	call tr2d(data_value(iad_x1),data_value(iad_y1),
     .            data_value(iad_x2),data_value(iad_y2),ns,
     .            o_alimet(1:1),o_alisig,ang,f1,f2,
     .            tr,xres,yres,mark,rms,ncl,stat,message)
	if(stat.ne.0)then
	  write(*,'(1x,a)')message
	  goto 9999
	endif

	do i=1,4
	  o_alitr(i) = tr(i)
	enddo
	o_alix0  = tr(5)
	o_aliy0  = tr(6)
	o_alif1  = f1
	o_alif2  = f2
	o_aliang = ang(1)

            write(*,*)
     1    '          f1        f2        ang      x0         y0'
            write(*,101)'fit: ',f1,f2,ang(1)*180/pi,tr(5),tr(6)
101         format(a,5(1x,f9.4))
            write(*,103)
     1      'erreur angle:',ang(2)*180/pi,' rms0:',rms0,' rms:',rms
103         format(3(a,f9.5))
            write(*,*)'  n  x-residue y-residue rejected'
            do 3000 k=1,ns
              write(*,102)k,xres(k),yres(k),mark(k)
102           format(1x,i2,2(1x,f9.6),6x,i1)
3000        continue
            write(*,*)ncl,'  clipping cycles'

	
	erreur_prog = 0.
9999	call reveil_inter(erreur)
	return
	end




        subroutine
     1  tr2d(x1,y1,x2,y2,nin,meth,fsig,ang,f1,f2,tr,xres,yres,mark,rms,
     2  ncl,stat,message)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c.identification:
c  subroutine trans		version 2.10	850812
c  k. banse			eso - garching
c  modif a. blecha              obsge           880822      
c    methode `l` ajoutee
c    sigma clipping
c    calcul des residues
c
c
c.purpose:
c  calculate the coefficients of a linear transformation between the
c  two point sets x2,y2 and x1,y1 of n points each.
c  the transformation is assumed to be of the form:
c	
c  x2(k) = tr(1)*x1(k) + tr(2)*y1(k) + tr(5)
c  y2(k) = tr(3)*x1(k) + tr(4)*y1(k) + tr(6)
c	

c.algorithm:
c  use least square fit to above equations
c
c  then compute parameters
c	
c  x2(k) =  f1*cos(ang)*x1(k) + f2*sin(ang)*y1(k) + tr(5)
c  y2(k) = -f1*sin(ang)*x1(k) + f2*cos(ang)*y1(k) + tr(6)
c
c.input/output:
c
c  input par:
c  ----------
c  x1:		r*4 array   2. set of points
c  y1
c  x2:		r*4 array   1. set of points
c  y2
c  nin:		i*4	    no. of points in x2 and x1
c  meth:	char.exp.   method to use,
c			    f = free paramaters
c                           e = f1 and f2 are equal
c			    u = f1 = f2 = 1.0
c			    s = shift only: f1 = f2 = 1.0, angle = 0.0
c	                    l = locked scale to value of f1: f2 = f1
c  fsig         r*4         factor of sigma clipping: reject cond.:
c                           fsig < 0 : residuals > -fsig*RMS
c                           fsig > 0 : residuals >   fsig
c                           fit is repeated until one of conditions
c                           is met:
c                           convergence: no new points rejected in two
c                                        consequitive clipping
c                           fit fails (not enough of points)
c                           max number (nclmax) of cycles is reached 
c
c  input-output or output only:
c  ---------------------
c  f1:		r*8		scaling factor in x-direction
c  f2:		r*8		scaling factor in y-direction
c
c  output only:
c  -----------
c  ang(2):	r*8		rotation angle in radians and its rms
c  tr:		r*8 array	array of coefficients
c				ang, f1, f2 are derived from tr ...
c  xres:        r*4 array 's    set of residuals x2(x1,y1)-x2
c  yres:                             and y2(x1,y1)-y2
c  mark:        car*1 array     set on 1 for points rejected during sigma
c  rms:                         rms of the transform in target units
c  stat:	i*4		return status
c
c-----------------------------------------------------------------------
c
c	
	character*(*)	meth,message
c	
	integer*4	i,nin,n,stat,mark(1)
c	
	real*4		x2(nin),y2(nin),x1(nin),y1(nin)
	real*4		xres(nin),yres(nin),fsig
c	
	real*8		tr(6),f1,f2,ang(1),ac,rms
	real*8		rn,x,y,xi,yi,xx,xy,yy,xxi,xyi,xiy,yyi
	real*8		a,b,c,d,dinv
	real*8		datan2,dcos,dsign,dsin,dsqrt
        data nclmax0/10/
	
	write(*,*)'x2 ',x2
	write(*,*)'y2 ',y2
	write(*,*)'x1 ',x1
	write(*,*)'y1 ',y1
	write(*,*)' '
c
c    set SIGMA clipping parameters
        ncl=0
        if (fsig.ne.0) then
          nclmax=nclmax0
        else
          nclmax=0
        endif
        do j=1,nin
          mark(j)=0
        enddo
        n=nin
c	
c  check input
1    	if ( ((meth(1:1).eq.'f').and.(n.lt.3)) .or.
     +	   ((meth(1:1).eq.'e'.or.meth(1:1).eq.'l').and.(n.lt.2)) .or.
     +	   ((meth(1:1).eq.'u').and.(n.lt.2)) .or.
     +	   ((meth(1:1).eq.'s').and.(n.lt.1)) ) then
	   message='not enough points in tables...'
	   stat = -1
	   return
	else
	   stat= 0
	endif
c	
c  init sums
	rn = 1.d0 / dble (float(n) )
	x = 0.d0
	y = 0.d0
	xi = 0.d0
	yi = 0.d0
	xx = 0.d0
	xy = 0.d0
	yy = 0.d0
	xxi = 0.d0
	xyi = 0.d0
	xiy = 0.d0
	yyi = 0.d0
c	
c  compute the sums
	do i=1,nin

          if (mark(i).eq.0) then
	   a = dble (x1(i))
	   b = dble (y1(i))
	   c = dble (x2(i))
	   d = dble (y2(i))
c	
	   x = x + a		
	   y = y + b
	   xx = xx + a*a
	   xy = xy + a*b
	   yy = yy + b*b
	   xi = xi + c
	   yi = yi + d
	   xxi = xxi + a*c
	   xyi = xyi + a*d
	   xiy = xiy + c*b
	   yyi = yyi + b*d
          endif

	enddo
c	
c  normalize coefficients -> all terms linear in x1 or y1 become 0
	xx = xx - x*x*rn
	xy = xy - x*y*rn
	yy = yy - y*y*rn
	xxi = xxi - x*xi*rn
	xyi = xyi - x*yi*rn
	xiy = xiy - xi*y*rn
	yyi = yyi - y*yi*rn
c	
c  branch according to method
	if (meth(1:1).eq.'e') goto 2000
	if (meth(1:1).eq.'u') goto 3000
	if (meth(1:1).eq.'s') goto 4000
	if (meth(1:1).eq.'l') goto 5000
c	
c  fall through to method = free	(all paramters are free)
	d = xx*yy - xy*xy
	if (dabs(d).lt.1.d-20) then
	   message='points not well chosen...'
	   stat = 1
	   return
	endif
c
	a = (yy*xxi - xy*xiy) / d
	b = (xx*xiy - xy*xxi) / d
	c = (yy*xyi - xy*yyi) / d
	d = (xx*yyi - xy*xyi) / d
c	
	tr(1) = a
	tr(2) = b
	tr(3) = c
	tr(4) = d
	tr(5) = (xi - a*x - b*y)*rn
	tr(6) = (yi - c*x - d*y)*rn
c	
cab     on prends la moyenne de deux determination d'angle
cab        et le rms comme indicateur d'erreur
cab
        ang1 = datan2(b,d)
	ang2 = datan2(-c,a)
        ang(1) = (ang1+ang2)/2
        ang(2) = sqrt((ang1*ang1+ang2*ang2)/2-ang(1)*ang(1))
	f1 = dsqrt(a*a + c*c)
	f2 = dsqrt(b*b + d*d)
	ac = dcos(ang(1))
	if (dabs(ac).gt.0.5d0) then
	   f1 = dsign( f1,tr(1)*ac)
	   f2 = dsign( f2,tr(4)*ac)
	else
	   ac = dsin(ang(1))
	   f1 = dsign( f1,-tr(3)*ac)
	   f2 = dsign( f2,tr(2)*ac)
	endif
	go to 6000
c	
c  method = equal	(f1 = f2)
2000	d = xx + yy
	if (dabs(d).lt.1.d-20) then
	   message='points not well chosen...'
	   stat = 1
	   return
	else
	   dinv = 1.d0 / d
	endif
c
	a = (xxi + yyi) * dinv
	b = (xiy - xyi) * dinv
c	
	tr(1) = a
	tr(2) = b
	tr(3) = -tr(2)
	tr(4) = tr(1)
	tr(5) = (xi - a*x - b*y)*rn
	tr(6) = (yi + b*x - a*y)*rn
c	
	ang(1) = datan2(b,a)
	f1 = dsqrt(a*a+b*b)
	f2 = f1
	go to 6000
c	
c  method = unit		(f1 = f2 = 1.0)
3000	a = datan2( (xiy-xyi),(xxi+yyi) )
	ang(1) = a
	tr(1) = dcos(ang(1))
	tr(2) = dsin(ang(1))
	tr(3) = -tr(2)
	tr(4) = tr(1)
	tr(5) = (xi - tr(1)*x - tr(2)*y)*rn
	tr(6) = (yi + tr(2)*x - tr(1)*y)*rn
c	
	f1 = 1.0d0
	f2 = 1.0d0
	go to 6000
c	
c  method = shift		(f1 = f2 = 1.0, angle = 0.0)
4000	ang(1) = 0.0d0
	tr(1) = 1.0d0
	tr(2) = 0.0d0
	tr(3) = 0.0d0
	tr(4) = 1.0d0
	tr(5) = (xi - x) * rn
	tr(6) = (yi - y) * rn
c	
	f1 = 1.0d0
	f2 = 1.0d0
c	
c  method = locked, echelle connue et identique X et Y = f1: f2 = f1)
5000	a = datan2( (xiy-xyi),(xxi+yyi) )
	ang(1) = a
	tr(1) = f1*dcos(ang(1))
	tr(2) = f1*dsin(ang(1))
	tr(3) = -tr(2)
	tr(4) = tr(1)
	tr(5) = (xi - tr(1)*x - tr(2)*y)*rn
	tr(6) = (yi + tr(2)*x - tr(1)*y)*rn
c	
	f2 = f1
	go to 6000
cab
cab     compute residuals
cab
6000    s2=0
        s20=0
        do j=1,nin
          xres(j)=tr(1)*x1(j)+tr(2)*y1(j)+tr(5)-x2(j)
          yres(j)=tr(3)*x1(j)+tr(4)*y1(j)+tr(6)-y2(j)
          d2=xres(j)**2+yres(j)**2
          s20=s20+d2
          if (mark(j).ne.1) then
            s2=s2+d2
          endif
        enddo
        rms=sqrt(s2/n)

        if (fsig.ne.0) then
c
c    nextsw - switch indiccating that something has been changed
c      from the previous itteration
          nextsw=0
          if (fsig.gt.0) then
            dsig=fsig**2
          else
            dsig=-fsig*s20/nin
          endif
          n    = 0

          do j=1,nin
            d2=xres(j)**2+yres(j)**2
            if (d2.gt.dsig) then
              if (mark(j).ne.1) nextsw=1
              mark(j)=1
            else
              if (mark(j).ne.0) nextsw=1
              mark(j)=0
              n=n+1
            endif
          enddo
          ncl=ncl+1
          if (nextsw.eq.0) return
          if (ncl.gt.nclmax) then
	    stat = 2
            message='Clipping does not converge'
            return
          endif
          go to 1
        endif
        return

c
	end