c**************************
c  Programme bmg - UNIX
c*************************************************************************
      subroutine bmg(nw,win,nm,iux,iuy,iug,gname,istop)
c                              iux inutilise
       character comm*80,fin*3,end*3
       character*22 gname
       dimension kz(512)
       integer seuili,seuilgi,seuilj,sgj,sgj1,sgj2,distor
       integer win(8),winp
c       integer displ
       integer*2 ima(1536,1536)
c                     attention images < 1536*1536
c---------------------
c                       voir SRECT
c      is=1536
c      js=1536
c     leps    ! intervalle max entre point probable (seuil) et gradient max
c      milanglei=60
c        anglei=milanglei*0.001
c      idangle=0.2*(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   ',ier)
      read(8,1)comm
      write(3,2)comm
        
        do n=1,8

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

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

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

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

       n2=nm-n1+1
       nc=(n1+n2)/2
       n3n4=0

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

30    sgj1=sgj  
      sgj2=sgj

       rewind(iuy)
       read(iuy)kz
       write(3,3)(kz(i),i=1,8)
       im=kz(2)
       jm=kz(3)
c
      write(3,*)' bmg: im,jm ',im,jm
          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,
     4           iux,iuy,iug,gname,
     5           ima,milgeo,istop)
           write(3,5)

       return
       end    

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,
     4                 iux,iuy,iug,gname,
     5                 ima,milgeo,istop)
      integer*2 lec(1536),listj(1536),ima(1536,1536),lecx(1536)
      dimension z(1536),zg(1536),x(1536),zmoy(1536),zmoyg(1536),
     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(1536,3),y(1536),
     3          vi(11,6,6),vj(11,6,6),flech(11,2),ili(11,2),
     4          xdes(7),ydes(7),avi(6,6),bvi(6,6),avj(6,6),bvj(6,6),
     5          kz(512),zd(1536),xbord(11,5,2),ybord(11,5,2),
     6          integ(2,    3)
c                   deb-fin
      integer si,sgi,sj,sgj1,sgj2,distor
c      integer displ
      character titre*45,gname*22,date*8,nomps*24

c          indices      xbord (X//i)
c                       ybord (y//j)
c                  (canal n)
c                 
c                                  n,5,1   n,5,2
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,4,1   n,4,2
      is=1536
      js=1536
      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


        call date_and_time(date)

      xa=0.1
      xb=0.6
      xc=0.7
      xd=0.9
      xcd=0.5*(xc+xd)

        do2 n=1,11
c                    10 bords detectees sur 4 cotes (+2 inutilises)
        do kg=1,2
          do ihb=1,5
          xbord(n,ihb,kg)=0.
          ybord(n,ihb,kg)=0.
          enddo
c                     6 points de reference sur 2 cotes
          do ihb=1,3
          xr(n,ihb,kg)=0.
          yr(n,ihb,kg)=0.
          enddo
        enddo
2       continue

      i2=im-i2m
      idangle=0.2*(i2-i1)
          write(3,*)'idangle',idangle
      ic=0.5*(i1+i2)+0.5
      centre=0.5*(i1+i2)
      icentre=ic  

      j2=jm-j2m
c      ia=ic-intvi
c      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,1536
        do ihmb=1,3
        zmi(j,ihmb)=0.
        enddo
        enddo
c

      rewind(iuy)
      read(iuy)kz
c---------------------          stockage de l'image entiere

       rewind(iux)
       read(iux)kz
       write(3,*)(kz(i),i=1,8)
c       im=kz(2)
c       jm=kz(3)

      iia=i1
      iib=i1+4
      iic=i2-4
      iid=i2

      do j=1,jm
      read(iux)(lecx(i),i=1,im)
      read(iuy)(lec(i),i=1,im)
       do i=1,im
       lec(i)=lec(i)-lecx(i)
       enddo

      if(j.eq.200)write(3,*)' bmg:lecx ',(lecx(i),i=1,100,10)
      if(j.eq.200)write(3,*)' bmg:lec  ',(lec(i),i=1,100,10)

c      call zero3(iia,iib,iic,iid,0,0,0,lec,1,im,1,0,
c     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
        ima(i,j)=lec(i)
        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

c3     continue
      write(3,*)' fond,zmax ',fond,zmax

      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',200.,20,'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
      if(iangle.eq.2)call pgsls(2)
      call pgline(jm,x,zd)
      call pgsls(1)
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',200.,20,'bcmts',50.,5)
      call pgwindow(x1,x2,z1,z2)
      jm1=jm-1
      call pgline(jm1,x,zg)
c
      zgseuil=sgj1*zgmax/100.  

      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',200.,20,'bcnts',200.,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
        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),1,0)
c                                            kg droite
        write(3,*)' droite',yc(n2,2,iangle)

        isig=+1
        icroiss=1
         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),2,0)
c                                            kg gauche
        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
          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
          ideb=max0(ideb,i1)
          ifin=min0(ifin,i2)
        isig=-1
        icroiss=0
      call maxgrad(z,ideb,ifin,isig,icroiss,zseuil,zgseuil,
c          -------
     1             yc(n,2,iangle),3,n)

        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
          ideb=max0(ideb,i1)
          ifin=min0(ifin,i2)

        isig=1
        icroiss=0
      call maxgrad(z,ideb,ifin,isig,icroiss,zseuil,zgseuil,
c          -------
     1             yc(n,1,iangle),4,n)

        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-----
      write(3,*)' End loop angle'
      type *,' End loop angle'

      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=0.5*(i1+i2)+0.5
      enddo
      enddo
c-------------------------------------------------------------------
c                           Bords en I    xbord(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 ---

      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 et (i1+i2)*0.5 

       dj=(0.5*(i1+i2)-ic)*anglei
c                                    voisin de 0
      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.

      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
c                                       xbord,ybord (n,4a5,kg)
260   continue 
c      write(3,*)' i1,icentre,i2 ',i1,icentre,i2
c                        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
c      if((zmoy(i+leps)-zseuili)*igsign.lt.0.)goto26
c      if((zmoy(i-leps)-zseuili)*igsign.gt.0.)goto26
        call SMAX(zmoyg,i,eps)
      piv=i+eps-0.5

      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)=piv
      y(2)=y(1)

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

      call pgline(2,x,y)
      goto27
26    continue
        Idessin=Idessin+1
27    continue
c
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
c      if((zmoy(i+leps)-zseuili)*igsign.gt.0.)goto28
c      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)
      piv=i+eps-0.5

      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)=piv
      y(2)=y(1)

         ybord(n,5,kg)=0.5*(x(1)+x(2))
         xbord(n,5,kg)=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                       -------
c
c                         Bords en J  xbord,ybord   (n,1a3,kg)
c                                         pour paraboles eventuelles

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

      do58 kg=1,2
c                                                 boucle kg
         call pgslw(3)
c      if(nkbon(n,kg).ne.1)goto58
c      if(xr(n,1,kg).eq.0..or.xr(n,3,kg).eq.0.)goto58

      xcc=0.5*(xbord(n,4,kg)+xbord(n,5,kg))
        icc=xcc
      xl=xbord(n,5,kg)-xbord(n,4,kg) 
      ili(n,kg)=lip*xl/100.

        do ihb=1,3
        integ(1,ihb)=icc+ili(n,kg)*(ihb-2)-intvi
        integ(2,ihb)=integ(1,ihb)+2*intvi
        enddo

      jdeb=yca(n,kg)-lj
      jfin=yca(n,kg)+lj
c
      jdeb=max0(jdeb,1)
      jfin=min0(jfin,jm)

51    format(' n='i2,'  kg='i2,' jdeb='i4,'   jfin='i4)
      write(3,51)n,kg,jdeb,jfin

c

52    do j=jdeb,jfin
        do ihb=1,3
           zmi(j,ihb)=0
           ip1=integ(1,ihb)
           ip2=integ(2,ihb)
              do ip=ip1,ip2
              yy=j+(ip-xcc)*anglei
              jp=yy
              dy=yy-jp
              zmi(j,ihb)=zmi(j,ihb)+ima(ip,jp)*(1.-dy)+
     1                              ima(ip,jp+1)*dy
              enddo
        enddo
      enddo
      do ihmb=1,3
      zmi(j,ihmb)=zmi(j,ihmb)/den
      enddo
c
      do54 ihb=1,3
      piv=0.
        do jp=jdeb+2,jfin-1
        zmoyg(jp)=zmi(jp,ihb)-zmi(jp-1,ihb)
        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)

      ybord(n,ihb,kg)=jmax+eps+
     1                (ili(n,kg)*(ihb-2)+xcc-icc)*anglei
      xbord(n,ihb,kg)=ili(n,kg)*(ihb-2)+xcc
c                                              flottants equidistants de xcc

53    format(' n,kg,ihb,  jmax,eps,ili,xcc '3i3,2x,i5,f7.2,i5,f6.1)
      write(3,53)n,kg,ihb,jmax,eps,ili(n,kg),xcc

        y(1)=integ(1,ihb)
        y(2)=integ(2,ihb)
c                                 dessin
             x(1)=ybord(n,ihb,kg)-intvi*anglei
             x(2)=ybord(n,ihb,kg)+intvi*anglei
535   format(' x1,x2,y1,y2 ',4f6.1)
      write(3,535)x(1),x(2),y(1),y(2)

      call pgline(2,x,y)
c          ------
54    continue
c              fin boucle ihb
58    continue
c                  fin boucle kg
60    continue
c                      fin boucle n
c-------------------------------------------------
      if(n3n4.ne.0)then
      call interpob(xbord,n3,n4)
      call interpob(ybord,n3,n4)
c          --------
      endif

611   format(/' xbord(n,ihb,kg) 11,5,2 ')
      write(3,611)
612   format(11f6.1)
      write(3,612)xbord
613   format(/' ybord(n,ihb,kg) 11,5 2 ')
      write(3,613)
      write(3,612)ybord

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

c                     on dispose de 10 points xbord,ybord centres
c                        ( xbord(2)=0.5*(xbord(1)+xbord(3)) )
c
c                     calcul de 6 points references xr,yr
c                     -----------------------------------
      call pgslw(1)
      do n=1,nm
c                                              calcul et dessin, cas lineaire
      if(distor.eq.0)then
      call interseclin(xbord,ybord,n,xr,yr)
c          -----------

      do63 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))

         do ihmb=1,3
         xx=yr(n,ihmb,kg)
         yy=xr(n,ihmb,kg)
         call pgpoint(1,xx,yy,5)
c             -------
         enddo
63      continue
      else
c                                              calcul et dessin, cas parabole
c                          on remplace xbord,ybord (1) et (3)
c                          par xbord,ybord (1) et (3) en (bas) et (haut)

      do66 kg=1,2
      x1=xbord(n,4,kg)
      x2=xbord(n,1,kg)
      x3=xbord(n,2,kg)
      x4=xbord(n,3,kg)
      x5=xbord(n,5,kg)
        rap=(x5-x1)/(x4-x2)
c      write(3,*)' n,kg,rap ',n,kg,rap

      y3=ybord(n,2,kg)
      y2=ybord(n,1,kg)
      y4=ybord(n,3,kg)
      dy=0.5*(y2+y4-2.*y3)*rap**2
      y1=y3+(y2-y4)*rap/2.+dy
      y5=y3+(y4-y2)*rap/2.+dy
c         write(3,*)' y1,y2,y3,y4,y5 ',y1,y2,y3,y4,y5

      xbord(n,1,kg)=x1
      xbord(n,3,kg)=x5

      ybord(n,1,kg)=y1
      ybord(n,3,kg)=y5
66    continue

      call interseclin(xbord,ybord,n,xr,yr)
c          -----------

      do kg=1,2
      xr(n,2,kg)=0.5*(xr(n,1,kg)+xr(n,3,kg))
      dxcc=xr(n,2,kg)-x3
      yr(n,2,kg)=ybord(n,2,kg)+dxcc*anglei
      enddo

      endif
c                       fin calcul xr,yr
      do kg=1,2
          flech(n,kg)=yr(n,2,kg)-0.5*(yr(n,1,kg)+yr(n,3,kg))
      enddo

c                          dessin contour canaux

      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)

      enddo
c                         fin boucle n
64    format(/'  XR(N,IHMB,KG)'/11f8.1/11f8.1/11f8.1//
     1                          11f8.1/11f8.1/11f8.1/)
      write(3,64)xr
65    format(/'  YR(N,IHMB,KG)'/11f8.1/11f8.1/11f8.1//
     1                          11f8.1/11f8.1/11f8.1/)
      write(3,65)yr
68    format(/'  Fleches(n,kg)'/11f8.1/11f8.1/)
      write(3,68)flech
      write(7,68)flech
c-------------------------------
70    format(3i8,a8,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')

78    format('  xr(n,ihmb,kg)  / yr(n,ihmb,kg)     rough values')
c      write(iug,70)kj,kh,kr,jan,mois,jour,gnom
      write(iug,70)kj,kh,kr,date,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

      write(iug,78)
      write(iug,72)xr
      write(iug,72)yr

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,64)xr
      write(3,65)yr
      write(iug,94)
      write(iug,72)xr
      write(iug,72)yr
c
95    continue

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
c      kgs=2
c      if(mark.eq.-1)kgs=1
c      khbs=khb+1
c      if(khb.eq.1)khbs=1
c

      rewind(iuy)  

110   continue
      write(iug,99)
      write(iug,72)ecart,angle1,angle2

120   continue
c                             dessin titre
        call date_and_time(date)
      write(titre(1:22),'(a22)')gname
      write(titre(23:27),'(a1)')'     '
      write(titre(28:35),'(a)')date
      write(titre(36:45),'(f10.1)')ecart
      call pgvport(xa,xb,.1,.9)
      call pgsch(0.7)
      call pgslw(2)
      call pglabel(' ',' ',titre)
c      close(unit=7)
           istop=0
           piv=float(milgeo)/1000.
           if(ecart.gt.piv)then
             type *,' Error too large '
             type *,' Check geo.ps'
             type *,' Modify the parameters or the accuracy threshold'
             istop=1
           endif
200   continue
      return
      end
c--------------------------
      subroutine SMAX(z,i,eps)
      dimension z(1)
      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(1536)
        do i=1,1536
        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,index,n)
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(1536),zg(1536)
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
        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
11    format(' Edge of channel not detected, index=',i1,
     1                                       ' n=',i2)
      if(zmax.eq.0.)then
        type 11,index,n
        write(3,11)index,n
        x=0.
        return
      endif

12    continue

      call smax(zg,imax,eps)
      x=imax+0.5+eps
        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(1536)

      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
        yy=j+(i-xc)*anglei
        jp=yy
        dy=yy-jp
        zp=zp+ima(i,jp)*(1.-dy)+ima(i,jp+1)*dy
        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(1536)

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.
        y1=j1+(i-xc)*anglei
        y2=j2+(i-xc)*anglei
        jp1=y1
        jp2=y2
        dy=y1-jp1
          do jp=jp1,jp2
          zp=zp+ima(i,jp)*(1.-dy)+ima(i,jp+1)*dy
          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,5,2),yb(11,5,2),xr(11,3,2),yr(11,3,2)

      call inters(xb,yb,n,1,1,3,1,4,1,4,2,xr,yr,1,1)
      call inters(xb,yb,n,1,1,3,1,5,1,5,2,xr,yr,3,1)
      call inters(xb,yb,n,1,2,3,2,4,1,4,2,xr,yr,1,2)
      call inters(xb,yb,n,1,2,3,2,5,1,5,2,xr,yr,3,2)

      return
      end
c--------------------------------
      subroutine inters(xb,yb,n,l1,k1,l2,k2,l3,k3,l4,k4,xr,yr,lr,kr)
c                ------
      dimension xb(11,5,2),yb(11,5,2),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,5,2)

      den=n4-n3
        do kg=1,2
        do ihmb=1,5
          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









