c**************************
c  Programme bmg - UNIX
c*************************************************************************
      subroutine bmg(nw,win,iux,iuy,iug,gname)
c                                  iux inutilise
       character comm*80,fin*3,end*3
       character*20 gname
       dimension kz(512)
       integer seuili,seuilgi,seuilj,sgj,sgj1,sgj2,distor
       integer win(8),winp
c       integer displ
       integer*2 ima(400,400)
c                     attention images < 400*400
c---------------------
c                       voir SRECT
c      is=400
c      js=400
c     leps    ! intervalle max entre point probable (seuil) et gradient max
c      milanglei=60
c        anglei=milanglei*0.001
c      idangle=0.1*(i2-i1)
c              demi-disrtance entre les 2 coupes utilisees pour estimer 
c                                                          l'angle/i   
c                  coupes a ic+/-idangle, integration sur +/-intvi
c     lip         pourcentage ecarts en I pour parabole haut/bas (bords J)
c          ex: 30
c      intercan=6  ! nombre de pixels entre canaux <---------- modifiable
c     icp       inutilise:
c         coupes //j autour de ic=(i1+i2)/2
c     lj        pour la determination precise des points en J,
c               recherche +/-lj  autour des valeurs approchees
c     intervi  intervalle integration en i pour detection bords j
c     intervj                            j                      i
c--------------------------------------
1      format(a)
2      format(1x,a)
3      format(12i6)
4      format(1x,12i6)
5      format(//)
6      format(3a4,3i8)
7      format(2x,3a4,3i8)

      rewind(8)
      call SNAME(8,'bmg   ')
      read(8,1)comm
      write(3,2)comm
        
        do n=1,8

      read(8,3) winp,milang,i1,i2m,icp,j1,j2m,lip,lj,
     1          intervi,intervj,intercan
      write(3,4)winp,milang,i1,i2m,icp,j1,j2m,lip,lj,
     1          intervi,intervj,intercan

      if(winp.eq.win(nw))goto20
        enddo
c-------
20    rewind(8)
      call sname(8,'bmg1  ')

      read(8,1)comm
      write(3,2)comm

       do n=1,8
      read(8,3)winp, seuili,seuilgi,seuilj,sgj,leps,nm,
     1         n1,n2,nc,n3n4,distor
      write(3,4)winp,seuili,seuilgi,seuilj,sgj,leps,nm,
     1         n1,n2,nc,n3n4,distor

      if(winp.eq.win(nw))goto30
        enddo

30      inver=0
      sgj1=sgj  
      sgj2=sgj

      khb=0
      mark=0
c              supprimer mark?
      kj=0
      kh=0
      kr=0      
      dia=0
      dja=0
      did=0
      djd=0
      jcor1=0
      jcor2=0

       rewind(iuy)
       read(iuy)kz
       write(3,3)(kz(i),i=1,8)
       im=kz(2)
       jm=kz(3)
c
          write(3,5)
      call srect(milang,las,i1,i2m,im,icp,j1,j2m,jm,lip,lj,
c          -----
     1           intervi,intervj,intercan,
     2           seuili,seuilgi,seuilj,sgj1,sgj2,leps,
     3           nm,n1,n2,nc,n3n4,distor,khb,inver,mark,kj,kh,kr,
     4           iux,iuy,iug,gname,
     5           dia,dja,did,djd,jcor1,jcor2,ima)
           write(3,5)
       type *,'avant stop'
c990    stop
       end    
c                     (rajouter la correction dark current par la suite)                


c----------------------
      subroutine SRECT(milang,las,i1,i2m,im,icp,j1,j2m,jm,lip,lj,
     1                 intvi,intvj,intercan,
     2                 si,sgi,sj,sgj1,sgj2,leps,nm,n1,n2,nc,n3n4,
     3                 distor,khb,inver,mark,
     4                 kj,kh,kr,iux,iuy,iug,gname,
     5                 dia,dja,did,djd,jcor1,jcor2,ima)
      integer*2 lec(1024),listj(1024),lecp(1024),ima(400,400)
      dimension z(1024),zg(1024),x(1024),zmoy(1024),zmoyg(1024),
     1          xr(11,3,2),yr(11,3,2),yc(11,2,2),cis(3,2),yca(11,2),
     2          ypar(11,3,2),nkbon(11,2),zmi(1024,3),y(1024),
     3          vi(11,6,6),vj(11,6,6),flech(11,2),xl2(11,2),xli(11,2),
     4          xdes(7),ydes(7),avi(6,6),bvi(6,6),avj(6,6),bvj(6,6),
     5          ind(12,11),kz(512),zd(1024),xbord(11,3,4),ybord(11,3,4)
      integer si,sgi,sj,sgj1,sgj2,distor,dia,dja,did,djd
c      integer displ
      character titre*43,gname*20,jdate*10,nomps*24

c          indices      xbord (X//i)
c                       ybord (y//j)
c                  (canal n)
c                 
c                                  n,1,4   n,2,4
c                                  --+-------+--
c                                  |           |
c                           n,3,1  +           +  n,3,2
c                                  |           |
c                           n,2,1  +           +  n,2,2
c                                  |           |
c                           n,1,1  +           +  n,1,2
c                                  |           |
c                                  --+-------+--
c   
c                                  n,1,3   n,2,3
      is=400
      js=400
      n3=n3n4/10
      n4=n3n4-10*n3
        if(n3.eq.0)n3=n1
        if(n4.eq.0)n4=n1

c     leps    ! intervalle max entre point probable (seuil) et gradient max

        anglei=milang*0.001

c      idangle=0.2*(i2-i1)
c                  coupes a ic+/-idangle, integration sur +/-intvi
c     lip         pourcentage ecarts en I pour parabole haut/bas (bords J)
c      intercan= nombre de pixels entre canaux
c     icp         inutilise

      do m=1,12
      do n=1,11
      ind(m,n)=0
      enddo
      enddo
c
        call idate(mois,jour,jan)
1     format(' date:  '3i8/)
c      write(3,1)jan,mois,jour


      xa=0.1
      xb=0.6
      xc=0.7
      xd=0.9
      xcd=0.5*(xc+xd)
        do2 n=1,11
        do2 kg=1,2
          yc(n,kg,1)=0.
          yc(n,kg,2)=0.
          nkbon(n,kg)=0
        do2 ihb=1,3
        xr(n,ihb,kg)=0.
        yr(n,ihb,kg)=0.
        ypar(n,ihb,kg)=0.
2       continue
      i2=im-i2m
      idangle=0.1*(i2-i1)
          write(3,*)'idangle',idangle
      ic=0.5*(i1+i2)+0.5
      centre=0.5*(i1+i2)
      icentre=ic  

      j2=jm-j2m
      ia=ic-intvi
      ib=ic+intvi
         iap=ic-intvi
         ibp=ic+intvi
c      write(3,*)' im,jm,i2,ic,j2,ia,ib ',im,jm,i2,ic,j2,ia,ib
c
        do j=1,1024
        do ihmb=1,3
        zmi(j,ihmb)=0.
        enddo
        enddo
c

      rewind(iuy)
      read(iuy)kz
c---------------------          stockage de l'image entiere
      iia=i1
      iib=i1+4
      iic=i2-4
      iid=i2
       if(inver.eq.1)then
        iia=im-iia+1
        iib=im-iib+1
        iic=im-iic+1
        iid=im-iid+1
       endif

      do j=1,jm
      read(iuy)(lec(i),i=1,im)

c        if(j.eq.200)write(3,*)(lec(i),i=1,286)
      call zero3(iia,iib,iic,iid,0,0,0,lec,1,im,1,0,
     1          az1,az2,az3,az4,az5,az6,az7,az8,az9)
c                          on retranche une fonction lineaire joignant
c                          les fonds en debut et fin de ligne
c        if(j.eq.200)write(3,*)' apres zero'
c        if(j.eq.200)write(3,*)(lec(i),i=1,286)
        do i=1,im
        ip=i
          if(inver.eq.1)ip=im-i+1 
        ima(i,j)=lec(ip)
        enddo
      enddo
c         write(3,*)'ima', ima(1,1),ima(im,jm)
c--------------

c      if(displ.eq.1)call pgbegin(0,'?',1,1)
c      if(displ.eq.2)call pgbegin(0,'/xwin',1,1)
c      if(displ.eq.3)then
c        nomps=gname//'.ps'
c        call system('rm '//nomps)
c        call pgbegin(0,nomps//'/ps',1,1)
c      endif

      call pgsch(0.6)
c
        fond=100000.
        zmax=0.
      type *,'jm :',jm

      do24 iangle=1,2    !!!  boucle iangle------------
        ia=iap+(2*iangle-3)*idangle
        ib=ibp+(2*iangle-3)*idangle

      call scanj(ima,is,js,ia,ib,j1,j2,anglei,z)
c          -----

      if(iangle.eq.1)then  !!  on definit zmax,zfond a la premiere coupe  --
        do j=j1,j2
        fond=amin1(z(j),fond)
        zmax=amax1(z(j),zmax)
        enddo
      endif

3     continue
c      write(3,*)' fond,zmax ',fond,zmax
c
c      if(iangle.eq.1)call pgadvance

      call pgvport(xa,xb,.70,.90)
        call pgslw(1)
        call pgsch(0.8)
      x1=1
      x2=jm
      z1=0
      z2=(zmax-fond)*1.
      call pgwindow(x1,x2,0.,100.)
      call pgbox('bcts',100.,10,'bcnts',20.,2)
      call pgwindow(x1,x2,z1,z2)

      xdes(1)=j1
      xdes(2)=j1
      ydes(1)=z1
      ydes(2)=z2
      call pgsls(2)
      call pgline(2,xdes,ydes)
      xdes(1)=j2
      xdes(2)=j2
      call pgline(2,xdes,ydes)
      call pgsls(1)
        do j=1,jm
        x(j)=j
        zd(j)=z(j)-fond
        enddo
      call pgline(jm,x,zd)

c
      zseuil=sj*(zmax-fond)/100.  !!  pour dessin
 
      x(2)=jm
      y(1)=zseuil
      y(2)=y(1)
      call pgsls(2)
      call pgline(2,x,y)
      call pgsls(1)
      x(2)=2
c
        if(iangle.eq.1)zgmax=0.
        do j=2,jm
        zg(j-1)=z(j)-z(j-1)
        piv=abs(zg(j-1))
        if(iangle.eq.1)zgmax=amax1(zgmax,piv)  !!  premiere coupe
        enddo

c        write(3,*)' zgmax',zgmax
      z2=zgmax
      z1=-z2

      call pgvport(xa,xb,.50,.70)
      call pgwindow(x1,x2,-100.,100.)
      call pgbox('abcnts',100.,10,'bcmts',50.,5)
      call pgwindow(x1,x2,z1,z2)
      jm1=jm-1
      call pgline(jm1,x,zg)
c
      zgseuil=sgj1*zgmax/100.  !!   pour dessin

      call pgsls(2)
      x(2)=jm
      y(1)=zgseuil
      y(2)=y(1)
      call pgline(2,x,y)
      y(1)=-y(1)
      y(2)=y(1)
      call pgline(2,x,y)
      call pgsls(1)
      x(2)=2
c
c                                      Bords approches en J   YC(N,KG)
      xb1=xb+0.02
      xb2=xb1+0.1
      piva=0.15
      pivb=piva+0.1
      call pgvport(xb1,xb2,piva,pivb)
      call pgwindow(0.,1.,0.,1.)
      call pgtext(0.,0.,'A-D')
      piva=piva+0.10
      pivb=piva+0.1
      call pgvport(xb1,xb2,piva,pivb)
      call pgwindow(0.,1.,0.,1.)
      call pgtext(0.,0.,'B-E')
      piva=piva+0.1
      pivb=piva+0.1
      call pgvport(xb1,xb2,piva,pivb)
      call pgwindow(0.,1.,0.,1.)
      call pgtext(0.,0.,'C-F')
c
      call pgvport(xa,xb,.10,.40)
      y1=1.
      y2=im
      call pgwindow(x1,x2,y1,y2)
      call pgbox('bcnts',100.,10,'bcnts',100.,10)
      call pglabel('J','I',' ')
      xdes(1)=1
      xdes(2)=jm
      ydes(1)=i1
      ydes(2)=i1
      call pgsls(2)
      call pgline(2,xdes,ydes)
      ydes(1)=i2
      ydes(2)=i2
      call pgline(2,xdes,ydes)
      call pgsls(1)
c
      zseuil=zseuil+fond
      zgseuil=zgseuil
c        write(3,*)' zseuil,zgseuil ',zseuil,zgseuil

c                                    bords extremes
        isig=-1
        icroiss=2
      call maxgrad(z,j1,j2,isig,icroiss,zseuil,zgseuil,
c          -------
     1             yc(n2,2,iangle))
c                                            kg droite
c        write(3,*)' droite',yc(n2,2,iangle)

        isig=+1
        icroiss=1
c         write(3,*)'****j1*j2*isig*icroiss*',j1,j2,isig,icroiss
      call maxgrad(z,j1,j2,isig,icroiss,zseuil,zgseuil,
c          -------
     1             yc(n1,1,iangle))
c                                            kg gauche
c        write(3,*)' gauche',yc(n1,1,iangle)

      can=(yc(n2,2,iangle)-yc(n1,1,iangle)+intercan)/float(n2-n1+1)
c         largeur moy canaux
c          write(3,*)' can ',can
c                                    bords KG=2   (droite canaux)
c                      
      do n=n1,n2-1
      iprobable=yc(n2,2,iangle)-(n2-n)*can
      ideb=iprobable-leps
      ifin=iprobable+leps
        isig=-1
        icroiss=0
      call maxgrad(z,ideb,ifin,isig,icroiss,zseuil,zgseuil,
c          -------
     1             yc(n,2,iangle))

        if(yc(n,2,iangle).eq.0.)then
        yc(n,2,iangle)=iprobable
        write(3,*)' attention point yc: iprobable=',iprobable
        endif
      enddo

c                                    bords kg=1   (gauche canaux)
      do n=n1+1,n2
      iprobable=yc(n1,1,iangle)+(n-n1)*can
      ideb=iprobable-leps
      ifin=iprobable+leps
        isig=1
        icroiss=0
      call maxgrad(z,ideb,ifin,isig,icroiss,zseuil,zgseuil,
c          -------
     1             yc(n,1,iangle))

        if(yc(n,1,iangle).eq.0.)then
          write(3,*)'  attention point yc: iprobable=',iprobable
          yc(n,1,iangle)=iprobable
        endif
      enddo  

      do n=n1,n2
        piv=yc(n,2,iangle)
        dpiv=(ib-ia)*0.5*anglei
        x(1)=piv-dpiv
        x(2)=piv+dpiv
        y(1)=ia
        y(2)=ib
        call pgline(2,x,y)
      enddo

      do n=n1,n2
        piv=yc(n,1,iangle)
        x(1)=piv-dpiv
        x(2)=piv+dpiv
        call pgline(2,x,y)
      enddo
 
20    format(/'  n1='i2,'  n2='i2)
c      write(3,20)n1,n2
22    format(/'  YC(N,KG,iangle)'/11f7.1/11f7.1/11f7.1/11f7.1)
c      write(3,22)yc

24    continue   !!! fin boucle iangle-----
c

      if(idangle.ne.0) call sangle(yc,n1,n2,idangle,angle)
c                           ------
      write(3,*)' angle ',angle
      anglei=angle

      do kg=1,2
      do n=n1,n2
        yca(n,kg)=0.5*(yc(n,kg,1)+yc(n,kg,2))
c                                              en ic=im*icp/100.
      enddo
      enddo
c-------------------------------------------------------------------
c                           Bords en I    XR(N,IHMB,KG)
      
        call pgslw(3)
        Idessin=0
c
      do40 n=n1,n2  !!!  boucle canaux -------------
           if(n.gt.n3.and.n.lt.n4)goto40

      do38 kg=1,2   !!!  boucle gauche-droite ---
c      if(mark.eq.1.and.kg.eq.1)goto38
c      if(mark.eq.-1.and.kg.eq.2)goto38
      if(n.eq.n2.and.kg.eq.2.and.idessin.ne.2)idessin=1
      if(yca(n,kg).eq.0.)goto38

      nkbon(n,kg)=1

c   decalage en j entre ic=im*icp/100. et (i1+i2)*0.5 
       ic=im*icp/100.
       dj=(0.5*(i1+i2)-ic)*anglei
      jdeb=yca(n,kg)+dj
      jfin=jdeb+intvj
        if(kg.eq.2)then
        jdeb=jdeb-intvj
        jfin=jfin-intvj
c                      rapportes a 0.5*(i1+i2)
        endif

      call scani(ima,is,js,i1,i2,jdeb,jfin,anglei,zmoy)
c          -----

        zseuili=fond+si*(zmax-fond)/100.
        zgseuili=sgi*zgmax/100.
        igsign=1
        if(sgi.lt.0)igsign=-1
        zgseuili=abs(zgseuili)
c          write(3,*)' zseuili,zgseuili ',zseuili,zgseuili

      do i=i1+1,i2
      zmoyg(i)=zmoy(i)-zmoy(i-1)
      enddo
      goto260
c
c                                      dessin coupe  en  I
250   Idessin=2
      call pgsls(1)
      call pgslw(1)
      x(1)=yca(n,kg)+dj+(i1-icentre)*anglei
      x(2)=yca(n,kg)+dj+(i2-1-icentre)*anglei
      y(1)=i1
      y(2)=i2-1
      call pgline(2,x,y)  !!! ligne non trouvee
        do i=i1,i2-1
        y(i)=i
        enddo
        zcmin=10000.
        zcmax=0.
        zcgmin=10000.
        zcgmax=-10000.
c       do i=i1,i2-1
c       zcmin=amin1(zmoy(i),zcmin)
c       zcmax=amax1(zmoy(i),zcmax)
c       zcgmin=amin1(zmoyg(i),zcgmin)
c       zcgmax=amax1(zmoyg(i),zcgmax)
c       enddo
      call pgvport(xcd,xd,.10,.40)

        piv=1.9*si
      call pgwindow(0.,piv,y1,y2)
       xsi=si
      call pgbox('bcmts',xsi,1,'bc',0.,0)

        piv=1.9*zseuili
      call pgwindow(0.,piv,y1,y2)
      call pgline(i2-1,zmoy,y)
c                      --
      call pgsls(2)
        x(1)=zseuili
        x(2)=x(1)
        y(2)=im
        call pgline(2,x,y)
      call pgsls(1)
c
      call pgvport(xc,xcd,.10,.40)
        piv=1.9*sgi
      call pgwindow(-piv,piv,y1,y2)
        xsgi=sgi
      call pgbox('bcmts',xsgi,1,'bc',0.,0)

        piv=1.9*zgseuili
      call pgwindow(-piv,piv,y1,y2)
      call pgsls(2)
        y(2)=im/2
        x(1)=zgseuili*igsign
        x(2)=x(1)
        call pgline(2,x,y)
        y(1)=im/2
        y(2)=im
        x(1)=-x(1)
        x(2)=x(1)
        call pgline(2,x,y)
      y(1)=1
      y(2)=2
      call pgsls(1)

      call pgline(i2-1,zmoyg,y)
c                      --
      call pgslw(3)
      call pgvport(xa,xb,.10,.40)
      call pgwindow(x1,x2,y1,y2)
      goto38
c
c                                     Detection bords en I
260   continue 
c      write(3,*)' i1,icentre,i2 ',i1,icentre,i2
                        if(khb.eq.2)goto27  !!! bas
      do26 i=i1+2,icentre-1
      if((zmoyg(i)-zmoyg(i-1))*igsign.lt.0.)goto26
      if((zmoyg(i)-zmoyg(i+1))*igsign.lt.0.)goto26
      if((zmoyg(i)-zgseuili*igsign)*igsign.lt.0.)goto26
      if((zmoy(i+leps)-zseuili)*igsign.lt.0.)goto26
      if((zmoy(i-leps)-zseuili)*igsign.gt.0.)goto26
        call SMAX(zmoyg,i,eps)
      piv=i+eps-0.5
      xr(n,1,kg)=piv
      jdeb=yca(n,kg)+dj
      jfin=jdeb+intvj
        if(kg.eq.2)then
        jdeb=jdeb-intvj
        jfin=jfin-intvj
        endif
      x(1)=jdeb+(piv-centre)*anglei
      x(2)=jfin+(piv-centre)*anglei
      y(1)=xr(n,1,kg)
      y(2)=y(1)

         ybord(n,kg,3)=0.5*(x(1)+x(2))
         xbord(n,kg,3)=y(1)

      call pgline(2,x,y)
      goto27
26    continue
        Idessin=Idessin+1
27    continue
c
                        if(khb.eq.1)goto29   !!! haut
      do28 ip=icentre+2,i2-1
      i=icentre+i2+1-ip
      if((zmoyg(i)-zmoyg(i-1))*igsign.gt.0.)goto28
      if((zmoyg(i)-zmoyg(i+1))*igsign.gt.0.)goto28
      if((zmoyg(i)+zgseuili*igsign)*igsign.gt.0.)goto28
      if((zmoy(i+leps)-zseuili)*igsign.gt.0.)goto28
      if((zmoy(i-leps)-zseuili)*igsign.lt.0.)goto28
270   format(' z, zg '4f8.0,2x4f8.0)
c      write(3,270)zmoy(i-leps),zmoy(i),zmoy(i+leps),zseuili,
c     1            zmoyg(i-1),zmoyg(i),zmoyg(i+1),zgseuili
        call SMAX(zmoyg,i,eps)
      xr(n,3,kg)=i+eps-0.5
        if(xr(n,1,kg).ne.0.)xr(n,2,kg)=0.5*(xr(n,1,kg)+xr(n,3,kg))
      piv=xr(n,3,kg)
      jdeb=yca(n,kg)+dj
      jfin=jdeb+intvj
        if(kg.eq.2)then
        jdeb=jdeb-intvj
        jfin=jfin-intvj
        endif
      x(1)=jdeb+(piv-centre)*anglei
      x(2)=jfin+(piv-centre)*anglei
      y(1)=xr(n,3,kg)
      y(2)=y(1)

         ybord(n,kg,4)=0.5*(x(1)+x(2))
         xbord(n,kg,4)=y(1)

      call pgline(2,x,y)
      goto29
28    continue
        if(Idessin.eq.0)Idessin=1
29    continue
c      if(mark.eq.1.or.mark.eq.-1)goto250
c
        if(Idessin.eq.1)goto250
c
38    continue !!!  fin boucle gche-droite ---
40    continue !!!  fin boucle canaux ----------------------------

      if(n3n4.ne.0)call interpo(xr,n3,n4)
c                       -------

45    format(/'  XR(N,IHMB,KG)'/11f8.1/11f8.1/11f8.1//
     1                          11f8.1/11f8.1/11f8.1/)
      write(3,45)xr

c      istop=1
c      if(istop.eq.1)goto200

c
c                         Bords en J    YPAR(N,IHMB,KG)   pour paraboles
c      if(mark.eq.2)then
c        do n=n1,n2
c        do kg=1,2
c          xli(n,kg)=1.
c        do ihmb=1,3
c          ypar(n,ihmb,kg)=yca(n,kg)+(xr(n,ihmb,kg)-ic)*angle
c        enddo
c        enddo
c        enddo
c        goto600
c      endif
c

      do60 n=n1,n2
                   if(n.gt.n3.and.n.lt.n4)goto60
      write(3,*)' ninter, n ',ninter,n
         den=2*intvi+1

      do58 kg=1,2
         call pgslw(3)
      if(nkbon(n,kg).ne.1)goto58
      if(xr(n,1,kg).eq.0..or.xr(n,3,kg).eq.0.)goto58
      xcc=0.5*(xr(n,1,kg)+xr(n,3,kg))
        icc=ic
      xl=xr(n,3,kg)-xr(n,1,kg)
      xl2(n,kg)=xl/2.
      xli(n,kg)=lip*xl/100.
        ip3=xcc-intvi
        ip4=xcc+intvi
        ip1=ip3-xli(n,kg)
        ip2=ip4-xli(n,kg)
        ip5=ip3+xli(n,kg)
        ip6=ip4+xli(n,kg)

      jdeb=yca(n,kg)-lj
      jfin=yca(n,kg)+lj
c
51    format(' n='i2,'  kg='i2,' jdeb='i4,'   jfin='i4)
c      write(3,51)n,kg,jdeb,jfin

c      goto200
c

52    do j=jdeb,jfin

      if(j.gt.jm)goto200
      if(j.lt.1)goto200

c        enddo
c      endif
c
        do ip=ip1,ip2
          jp=j+(ip-icc)*angle
        zmi(j,1)=zmi(j,1)+ima(ip,jp)
        enddo
        do ip=ip3,ip4
          jp=j+(ip-icc)*angle
        zmi(j,2)=zmi(j,2)+ima(ip,jp)
        enddo
        do ip=ip5,ip6
          jp=j+(ip-icc)*angle
        zmi(j,3)=zmi(j,3)+ima(ip,jp)
        enddo
      enddo
      do ihmb=1,3
      zmi(j,ihmb)=zmi(j,ihmb)/den
      enddo
c
      do54 ihmb=1,3
      piv=0.
        do jp=jdeb+2,jfin-1
        zmoyg(jp)=zmi(jp,ihmb)-zmi(jp-1,ihmb)
        if(kg.eq.1)piv=amax1(piv,zmoyg(jp))
        if(kg.eq.2)piv=amin1(piv ,zmoyg(jp))
        if(zmoyg(jp).eq.piv)jmax=jp
        enddo
      if(abs(zmoyg(jmax)).lt.zgseuil)jmax=0
      if(jmax.eq.(jdeb+2).or.jmax.eq.(jfin-1))jmax=0
      eps=0.
      if(jmax.eq.0)goto54
      call SMAX(zmoyg,jmax,eps)

      ypar(n,ihmb,kg)=jmax+eps-0.5+(xr(n,ihmb,kg)-ic)*angle
     

53    format(' n,kg,ihmb  :  jmax,eps '3i3,2x,i5,f7.2)
c      write(3,53)n,kg,ihmb,jmax,eps
        if(ihmb.eq.1)then
        y(1)=ip1
        y(2)=ip2
             x(1)=yca(n,kg)+(ip1-icc)*angle
             x(2)=yca(n,kg)+(ip2-icc)*angle
        endif
        if(ihmb.eq.2)then
        y(1)=ip3
        y(2)=ip4
             x(1)=yca(n,kg)+(ip3-icc)*angle
             x(2)=yca(n,kg)+(ip4-icc)*angle
        endif
        if(ihmb.eq.3)then
        y(1)=ip5
        y(2)=ip6
             x(1)=yca(n,kg)+(ip5-icc)*angle
             x(2)=yca(n,kg)+(ip6-icc)*angle
        endif

        ybord(n,ihmb,kg)=0.5*(x(1)+x(2))
        xbord(n,ihmb,kg)=0.5*(y(1)+y(2))

      call pgline(2,x,y)
54    continue
58    continue
60    continue

      if(n3n4.ne.0)then
      call interpo(ypar,n3,n4)
      call interpob(xbord,n3,n4)
      call interpob(ybord,n3,n4)
c          --------
      endif

61    format(/'   YPAR(N,IHMB,KG)'/11f8.1/11f8.1/11f8.1//
     1                             11f8.1/11f8.1/11f8.1/)
      write(3,61)ypar

611   format(/' xbord, puis ybord (n,kg,ihmb)')
      write(3,611)
612   format(11f6.1)
      write(3,612)xbord,ybord

c
c                         cisaillement  cis(ihmb,kg)
c
      den=yca(nc,2)-yca(nc,1)-intvj
      if(den.eq.0.)goto620
        do ihmb=1,3,2
        alp=(xr(nc,ihmb,2)-xr(nc,ihmb,1))/den
        cis(ihmb,1)=-alp*intvj/2.
        cis(ihmb,2)=alp*intvj/2.
        enddo
      cis(2,1)=0.5*(cis(1,1)+cis(3,1))
      cis(2,2)=0.5*(cis(1,2)+cis(3,2))
62    format(' Cisaillements  ' 6f6.1/)
620   continue
      write(3,62)cis
c
c                         3 points references
      call pgslw(1)
      do n=1,nm
      do63 kg=1,2

      if(distor.eq.0)then
      call interseclin(xbord,ybord,n,xr,yr)
c          -----------
         do ihmb=1,3
         xx=yr(n,ihmb,kg)
         yy=xr(n,ihmb,kg)
         call pgpoint(1,xx,yy,5)
         enddo
      goto63
      endif

          flech(n,kg)=0.
        yy1=ypar(n,1,kg)
        yy2=ypar(n,2,kg)
        yy3=ypar(n,3,kg)
      if(distor.eq.0)yy2=0.5*(yy1+yy3)
      if(yy1*yy2*yy3.eq.0.)goto63
      cp=yy2
      bp=(yy3-yy1)/(2.*xli(n,kg))
      ap=(yy1+yy3-2.*yy2)/(2.*xli(n,kg)*xli(n,kg))
      yr(n,2,kg)=yy2               
      yr(n,1,kg)=cp-xl2(n,kg)*(bp-xl2(n,kg)*ap)
      yr(n,3,kg)=cp+xl2(n,kg)*(bp+xl2(n,kg)*ap)
        flech(n,kg)=yr(n,2,kg)-0.5*(yr(n,1,kg)+yr(n,3,kg))
c
        do ihmb=1,3
          xxr=xr(n,ihmb,kg)
          if(xxr.ne.0.)xr(n,ihmb,kg)=xxr+cis(ihmb,kg)
        xx=yr(n,ihmb,kg)
        yy=xr(n,ihmb,kg)
        call pgpoint(1,xx,yy,5)
        enddo
63    continue
c
c      if(n.eq.nc)then
      xdes(1)=yr(n,1,1)
      ydes(1)=xr(n,1,1)
      xdes(2)=yr(n,2,1)
      ydes(2)=xr(n,2,1)
      xdes(3)=yr(n,3,1)
      ydes(3)=xr(n,3,1)
      xdes(4)=yr(n,3,2)
      ydes(4)=xr(n,3,2)
      xdes(5)=yr(n,2,2)
      ydes(5)=xr(n,2,2)
      xdes(6)=yr(n,1,2)
      ydes(6)=xr(n,1,2)
      xdes(7)=xdes(1)
      ydes(7)=ydes(1)
      call pgsls(4)
      call pgline(7,xdes,ydes)
      call pgsls(1)
c      endif
      enddo
c
      write(3,45)xr
64    format(/'  YR(N,IHMB,KG)'/11f8.1/11f8.1/11f8.1//
     1                          11f8.1/11f8.1/11f8.1/)
      write(3,64)yr
66    format(/'  Fleches(n,kg)'/11f8.1/11f8.1/)
      write(3,66)flech
c
70    format(3i8,3i4,20xa15)
72    format(11f7.1)
74    format(12i6)
75    format(' displ   las    i1   i2m   icp    j1   j2m   lip    lj',
     1       ' intvi intvj')
76    format('    si   sgi    sj  sgj1   sgj2 leps    nm    n1    n2',
     1       '    nc distor')
77    format(' distor  khb inver  mark   DiA   DjA   DiD   DjD')
78    format('  xr(n,ihmb,kg)  / yr(n,ihmb,kg)     rough values')
      write(iug,70)kj,kh,kr,jan,mois,jour,gnom
      write(iug,75)
      write(iug,74)displ,las,i1,i2m,icp,j1,j2m,lip,lj,intvi,intvj
      write(iug,76)
      write(iug,74)si,sgi,sj,sgj1,sgj2,leps,nm,n1,n2,nc,distor
c      write(iug,77)
c      write(iug,74)distor,khb,inver,mark,dia,dja,did,djd,jcor1,jcor2
      write(iug,78)
      write(iug,72)xr
      write(iug,72)yr
      if(mark.ne.0)goto940
c
c                          vecteurs  en  fonction de N
      do n=1,nm
      do ihmb1=1,3
      do kg1=1,2
      ind1=ihmb1+3*(kg1-1)
      do ihmb2=1,3
      do kg2=1,2
      ind2=ihmb2+3*(kg2-1)
        vi(n,ind1,ind2)=0.
        u=xr(n,ihmb1,kg1)
        v=xr(n,ihmb2,kg2)
        if(u*v.ne.0.)vi(n,ind1,ind2)=v-u
        vj(n,ind1,ind2)=0.
        u=yr(n,ihmb1,kg1)
        v=yr(n,ihmb2,kg2)
        if(u*v.ne.0.)vj(n,ind1,ind2)=v-u
      enddo
      enddo
      enddo
      enddo
      enddo
c
80    format(' Vecteurs  X       ABCDEF * ABCDEF  * n '/)
81    format(' Vecteurs  Y       ABCDEF * ABCDEF  * n '/)
82    format(2x11f7.1/2x11f7.1/2x11f7.1/2x11f7.1/2x11f7.1/2x11f7.1//)
c      write(3,80)
c      write(3,82)vi
c      write(3,81)
c      write(3,82)vj
c
      call pgvport(xc,xcd,0.5,0.9)
      call pglabel(' ',' ','I')
      call pgvport(xcd,xd,0.5,0.9)
      call pglabel(' ',' ','J')
c
      xp1=xc
      xp2=xd
      yp1=0.50
      h=(.90-.50)/7.
      yp1=yp1-h
c
      call slin(vi,nm,avi,bvi,ecarti)
      call slin(vj,nm,avj,bvj,ecartj)
      ecart=amax1(ecarti,ecartj)
c
      angle1=0.
      piv=avj(1,4)+nc*bvj(1,4)
      if(piv.ne.0.)
     1 angle1=10000.*(avi(1,4)+nc*bvi(1,4))/piv
      angle2=0.
      piv=avj(3,6)+nc*bvj(3,6)
      if(piv.ne.0.)
     1 angle2=10000.*(avi(3,6)+nc*bvi(3,6))/piv
c
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,1,4,'AD',avi,bvi,avj,bvj)
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,2,5,'BE',avi,bvi,avj,bvj)
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,3,6,'CF',avi,bvi,avj,bvj)
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,1,3,'AC',avi,bvi,avj,bvj)
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,4,6,'DF',avi,bvi,avj,bvj)
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,1,2,'AB',avi,bvi,avj,bvj)
      call sgraph(xp1,xp2,yp1,h,vi,vj,nm,4,5,'DE',avi,bvi,avj,bvj)
c     
c                                             calcul points manquants 
      call pgvport(xa,xb,.10,.40)
      call pgwindow(x1,x2,y1,y2)
      do84 n=1,nm
      do84 ihmb=1,3
      do84 kg=1,2
      xpiv=xr(n,ihmb,kg)
      ypiv=yr(n,ihmb,kg)
      if(xpiv*ypiv.ne.0.)goto84
      kg2=3-kg
      xpiv2=xr(n,ihmb,kg2)
      ypiv2=yr(n,ihmb,kg2)
      if(xpiv2*ypiv2.eq.0.)goto84
      ind1=ihmb+3*(kg-1)
      ind2=ihmb+3*(kg2-1)
      vvi=avi(ind1,ind2)+n*bvi(ind1,ind2)
      vvj=avj(ind1,ind2)+n*bvj(ind1,ind2)
      xpiv=xpiv2-vvi
      ypiv=ypiv2-vvj
      xr(n,ihmb,kg)=xpiv
      yr(n,ihmb,kg)=ypiv
      call pgpoint(1,ypiv,xpiv,4)
84    continue
c
94    format('  xr / yr  with complementary points')
940   write(3,94)
      write(3,45)xr
      write(3,64)yr
      write(iug,94)
      write(iug,72)xr
      write(iug,72)yr
c
95    continue
96    format(' extreme marks  y,x(n1,1(or3),2) / y,x(n2,1(or3),1)')
196   format(' exreme marks 1st channel  ',
     1       'y,x(1,1,2(or1)) / y,x(1,3,2(or1))')
97    format(' arrows')
98    format(/)
99    format(' discrep.,angle1,angle2')
      write(3,98)
      write(3,99)
      write(3,72)ecart,angle1,angle2
      write(iug,97)
      write(iug,72)flech
c
      kgs=2
      if(mark.eq.-1)kgs=1
      khbs=khb+1
      if(khb.eq.1)khbs=1
c
      if(mark.eq.1.or.mark.eq.-1)then
c      write(3,196)
c      write(3,98)
c      write(3,72)yr(1,1,kgs),xr(1,1,kgs)
c      write(3,72)yr(1,3,kgs),xr(1,3,kgs)
      goto110
      endif
c
c      write(3,98)

c
      rewind(iuy)  
      write(3,96)
      write(3,72)yr(n1,khbs,2),xr(n1,khbs,2)
      write(3,72)yr(n2,khbs,1),xr(n2,khbs,1)
c
110   continue
      write(iug,99)
      write(iug,72)ecart,angle1,angle2
c
      if(mark.eq.1.or.mark.eq.-1)then
      write(iug,196)
      write(iug,72)yr(1,1,kgs),xr(1,1,kgs)
      write(iug,72)yr(1,3,kgs),xr(1,3,kgs)
      goto120
      endif
c
      write(iug,96)
      write(iug,72)yr(n1,khbs,2),xr(n1,khbs,2)
      write(iug,72)yr(n2,khbs,1),xr(n2,khbs,1)
c
112   format(' jA jD iA iD')
      write(iug,112)
      do n=1,11
      ind(1,n)=dja+yr(n,1,1)+0.5
      ind(4,n)=djd+yr(n,1,2)+0.5
      ind(7,n)=dia+xr(n,1,1)+0.5
      ind(10,n)=did+xr(n,1,2)+0.5
      enddo
      write(iug,74)ind
c
      write(3,112)
      write(3,74)ind
120   continue
        call date(jdate)
      write(titre(1:20),'(a20)')gname
      write(titre(21:25),'(a1)')'     '
      write(titre(26:35),'(a)')jdate
      write(titre(36:43),'(f8.1)')ecart
      call pgvport(xa,xb,.1,.9)
      call pgsch(0.7)
      call pgslw(2)
      call pglabel(' ',' ',titre)
c      close(unit=7)
200   continue
      return
      end
c--------------------------
      subroutine SMAX(z,i,eps)
      dimension z(4000)
      eps=0.
      b=z(i+1)-z(i-1)
      a=z(i+1)+z(i-1)-2.*z(i)
        if(a.eq.0.)return
      eps=-b/(2.*a)
      return
      end
c---------------------------
      subroutine sgraph(xp1,xp2,yp1,h,vi,vj,nm,ind1,ind2,tit,
     1                  avi,bvi,avj,bvj)
      dimension vi(11,6,6),vj(11,6,6),v(11),avi(6,6),bvi(6,6),
     1                                      avj(6,6),bvj(6,6),
     2          xdes(2),ydes(2)
      character tit*2
1     format(/' vecteurs  ind1='i1,'  ind2='i1,2x,11f7.1)
      yp1=yp1+h
      yp2=yp1+h
c
      do100 ij=1,2
      u1=xp1
      u2=xp2
      if(ij.eq.1)u2=0.5*(u1+u2)
      if(ij.eq.2)u1=0.5*(u1+u2)
      vmin=10000.
      vmax=-10000.
c
        do50 n=1,nm
        if(ij.eq.1)v(n)=vi(n,ind1,ind2)
        if(ij.eq.2)v(n)=vj(n,ind1,ind2)
        if(v(n).eq.0.)goto50
        vmin=amin1(vmin,v(n))
        vmax=amax1(vmax,v(n))
50      continue
c      write(3,1)ind1,ind2,v
c
      if(abs(vmax-vmin).gt.15.)goto100
      vmin=vmin-2
      vmax=vmax+2
      call pgvport(u1,u2,yp1,yp2)
      xnm1=nm+1
      call pgwindow(0.,xnm1,vmin,vmax)
      call pgbox('bcs',si,1,'bcs',10.,10)
      test=vmax-vmin
      if(test.lt.10)call pgbox(' ',1.,0,'bcs',10.,10)
c
        do60 n=1,nm
        if(v(n).eq.0.)goto60
        xx=n
        call pgpoint(1,xx,v(n),5)
60      continue
c
      xdes(1)=1.
      xdes(2)=nm
      if(ij.eq.1)then
        ai=avi(ind1,ind2)
        bi=bvi(ind1,ind2)
        ydes(1)=ai+bi
        ydes(2)=ai+nm*bi
      endif
      if(ij.eq.2)then
        aj=avj(ind1,ind2)
        bj=bvj(ind1,ind2)
        ydes(1)=aj+bj
        ydes(2)=aj+nm*bj
      endif
      call pgline(2,xdes,ydes)
100   continue
c
      u1=xp2
      u2=xp2+0.05
      call pgvport(u1,u2,yp1,yp2)
      call pgwindow(0.,1.,0.,1.)    
      call pgtext(0.4,0.4,tit)
      return
      end
c------------------------
      subroutine SLIN(v,nm,av,bv,ecart)
      dimension v(11,6,6),av(6,6),bv(6,6),
     1          x(11),y(11),p(11),coef(10)
      ecart=0.
      do4 i1=1,6
      do4 i2=1,6
      av(i1,i2)=0.
      bv(i1,i2)=0.
c
        nd=0
        do2 n=1,nm
        x(n)=n
        y(n)=v(n,i1,i2)
        p(n)=0.
        if(y(n).eq.0.)goto2
        p(n)=1.
        nd=nd+1
2       continue
c
      if(nd.lt.3)goto4
      call DPMCAR(x,y,p,nm,2,coef)
      av(i1,i2)=coef(1)
      bv(i1,i2)=coef(2)
c
      do3 n=1,nm
      yy=v(n,i1,i2)
      if(yy.eq.0.)goto3
      piv=av(i1,i2)+n*bv(i1,i2)-yy
      piv=abs(piv)
      ecart=amax1(ecart,piv)
3     continue
c
4     continue
      return
      end
c--------------------
      subroutine LISTEJ(listj)
      integer*2 listj(1024)
        do i=1,1024
        listj(i)=listj(i)/100
        enddo
1     format(i6,2x16i4)
2     format('     j')
c      write(3,2)
      do n=1,64
      i1=1+(n-1)*16
      i2=i1+15
c      write(3,1)i1,(listj(i),i=i1,i2)
      enddo
      return
      end
c---------------------
      subroutine sangle(yc,n1,n2,id,a)
c                ------
      dimension yc(11,2,   2)
c                     g/d  iangle
c      write(3,*)' n1,n2,id ',n1,n2,id

         a=0.
      do n=n1,n2
      a=a+yc(n,1,2)+yc(n,2,2)-yc(n,1,1)-yc(n,2,1)
      enddo

      a=a/(4.*(n2-n1+1)*float(id))
      return
      end
c-----------------------
      subroutine maxgrad(z,i1,i2,isig,icroiss,zseuil,zgseuil,x)
c                -------
c                          on cherche entre i1 et i2
c                               le point x de gradient max  de signe isig
c           icroiss=0   max absolu sur l'intervalle
c                   1   premier max dans le sens des i croissants
c                   2                                  decroissants
c           on ne retient que les points pour lesquels
c                                                     z>zseuil
c                                                 abs(zg)>zgseuil
      dimension z(1024),zg(1024)
c      write(3,*)' '
c      write(3,*)' isig,icroiss ',isig,icroiss
c      write(3,*)' i1,i2, zseuil,zgseuil ',i1,i2,zseuil,zgseuil

      do i=i1,i2-1
      zg(i)=z(i+1)-z(i)
      enddo
        zg(i2)=zg(i2-1)
c      do i=i1,i2-1
c       write(3,*)' i,z,zg', i,z(i),zg(i)
c      enddo

        zmax=0.

      do10 ip=i1+1,i2-1
        i=ip
        if(icroiss.eq.2)i=i2-ip+i1
c                              pour icroiss=0, on garde ip
        
      if(z(i).lt.zseuil)goto10
        piv2=isig*zg(i)
      if(piv2.lt.zgseuil)goto10
      if(piv2.lt.zgseuil)goto10
        piv1=isig*zg(i-1)
        piv3=isig*zg(i+1)
      if(piv2.lt.piv1.or.piv2.lt.piv3)goto10
c        write(3,*),' z, zg(i-1/i/i+1) ',z(i),zg(i-1),zg(i),zg(i+1)

      if(zmax.le.piv2)imax=i
      zmax=amax1(zmax,piv2)
      if(icroiss.ne.0)goto12
10    continue

      if(zmax.eq.0.)then
        type *,' pas de max dans maxgrad '
        x=0.
        return
      endif

12    continue

      call smax(zg,imax,eps)
      x=imax+0.5+eps
c        write(3,*)' x ',x
      return
      end
c----------------------
      subroutine scanj(ima,is,js,i1,i2,j1,j2,anglei,z)
c                -----
      integer*2 ima(is,js)
      dimension z(1024)

      den=i2-i1+1
      xc=0.5*(i1+i2)
c                         integ entre i1 et i2
      do j=j1,j2
      zp=0.
        do i=i1,i2
        jp=j+(i-xc)*anglei
        zp=zp+ima(i,jp)
        enddo
      
      z(j)=zp/den
      enddo
      return
      end
c------------------------
      subroutine scani(ima,is,js,i1,i2,j1,j2,anglei,z)
c                -----
      integer*2 ima(is,js)
      dimension z(1024)

c      write(3,*)' i1,i2,j1,j2 ',i1,i2,j1,j2

      den=j2-j1+1
      xc=0.5*(i1+i2)
c                         integ entre j1 et j2 autour de xc
      do i=i1,i2
      zp=0.
        jp1=j1+(i-xc)*anglei+0.5
        jp2=j2+(i-xc)*anglei+0.5
          do jp=jp1,jp2
          zp=zp+ima(i,jp)
          enddo

      z(i)=zp/den
      enddo
      return
      end
c-------------------------------
      subroutine interseclin(xb,yb,n,xr,yr)
c                -----------
c recherche des sommets xr,yr

      dimension xb(11,3,4),yb(11,3,4),xr(11,3,2),yr(11,3,2)

      call inters(xb,yb,n,1,1,3,1,1,3,2,3,xr,yr,1,1)
      call inters(xb,yb,n,1,1,3,1,1,4,2,4,xr,yr,3,1)
      call inters(xb,yb,n,1,2,3,2,1,3,2,3,xr,yr,1,2)
      call inters(xb,yb,n,1,2,3,2,1,4,2,4,xr,yr,3,2)
        do kg=1,2
        xr(n,2,kg)=0.5*(xr(n,1,kg)+xr(n,3,kg))
        yr(n,2,kg)=0.5*(yr(n,1,kg)+yr(n,3,kg))
        enddo
      return
      end
c--------------------------------
      subroutine inters(xb,yb,n,l1,k1,l2,k2,l3,k3,l4,k4,xr,yr,lr,kr)
c                ------
      dimension xb(11,3,4),yb(11,3,4),xr(11,3,2),yr(11,3,2)

      x1=xb(n,l1,k1)
      x2=xb(n,l2,k2)
      x3=xb(n,l3,k3)
      x4=xb(n,l4,k4)
c         write(3,*)' n,x1,x2,x3,x4 ',n,x1,x2,x3,x4

      y1=yb(n,l1,k1)
      y2=yb(n,l2,k2)
      y3=yb(n,l3,k3)
      y4=yb(n,l4,k4)
c         write(3,*)' n,y1,y2,y3,y4 ',n,y1,y2,y3,y4

        ayx=(y2-y1)/(x2-x1)
        byx=(y4-y3)/(x4-x3)
        axy=1./ayx
        bxy=1./byx

      xx=(x3*byx-x1*ayx+y1-y3)/(byx-ayx)
      yy=(y3*bxy-y1*axy+x1-x3)/(bxy-axy)

      xr(n,lr,kr)=xx
      yr(n,lr,kr)=yy

      return
      end
c-------------------------------------
      subroutine interpo(z,n3,n4)
c                -------
      dimension z(11,3,2)

      den=n4-n3
        do kg=1,2
        do ihmb=1,3
          do n=n3+1,n4-1
          x=float(n-n3)/den
          z(n,ihmb,kg)=(1.-x)*z(n3,ihmb,kg)+x*z(n4,ihmb,kg)
          enddo
        enddo
        enddo
      return
      end
  
c-------------------------------------
      subroutine interpob(z,n3,n4)
c                --------
      dimension z(11,3,4)

      den=n4-n3
        do kg=1,4
        do ihmb=1,3
          do n=n3+1,n4-1
          x=float(n-n3)/den
          z(n,ihmb,kg)=(1.-x)*z(n3,ihmb,kg)+x*z(n4,ihmb,kg)
          enddo
        enddo
        enddo
      return
      end









