c program qui trace les donnees sortant de SEXTRACTOR et ayant etant
c selectionnees comme target.

c* variables de lecture du fichier comportant les donnees*/
      integer i
      real xs(700),ys(700),as(700),bs(700),thetas(700)
c      character file_data
      
c* variables pour smongo */

      character smongo_data*80
      integer iplot
      character printer*20,title*100



c* calcul des ellipses
      integer i_ang, dimen_x1, dimen_x2, nb_tot
      real ang(550), x(550), y(550), div(550)
      real x1(550), x2(550), y1(550), y2(550)
      real a(710), b(710)
      real rot(550)
      real ang_0
      real PI
      real ellip
      character header*80

      PI = 3.1415927
      
      write(*,*)"entrer le nom du fichier smongo:"
      read(*,*)smongo_data
      open(unit=1,file=smongo_data)
      read(1,*)nbe
      do i=1, nbe
         read(1,*) xs(i), ys(i), as(i), bs(i), thetas(i)
      enddo


      write(*,600)
 600  format(2x,"Option graphique:  GTERM    (1)",/, 
     .        2x," LASER    (2)",/, 
     .        2x,"PSFILE   (3):", $)
      read(*,*) iplot

      if(iplot.eq.1) printer = 'xterm'
      if(iplot.eq.2) printer = 'postland ldouble2'
      if(iplot.eq.3) then
         printer = 'postportfile gr.ps'
      endif

      write(*,*)"Enter the title of the graph:"
      read(*,*)title

      call sm_device(printer)
      call sm_graphics
      call sm_erase
      call sm_defvar("TeX_strings","1")
      call sm_expand(1.001)
      call sm_window(1,1,1,1,1,1)
      call sm_location(6000,29000,10000,26200)
      call sm_limits(1.,1024.,1.,1024.)
      call sm_box(1,2,0,0)
      call sm_xlabel("Pixels")
      call sm_ylabel("Pixels")

      call sm_relocate(512.,1035.) 
      write(header,102) title
 102  format(a100)
      call sm_putlabel(8,header)

    
      do n=1,nbe
         thetas(n)= -thetas(n)
         as(n)=as(n)*10.
         bs(n)=bs(n)*10.
c pour mettre l'angle en radian car thetas est en degre
         rot(n) = thetas(n)*PI/180.

         dimen_x1=0
         dimen_x2=0

         ang_0=(-PI+0.1)/2.
c         ang_0 = -PI/2.
         do i_ang=1,50
            ang(i_ang)=0.
         enddo

         do i_ang=1,50
            ang(i_ang) = ang_0 + (i_ang-1) * PI/50.0
c            write(12,*)i_ang, ang(i_ang), ang_0, (i_ang-1)*PI/50.

            x(i_ang) = as(n)*cos(ang(i_ang))
            y(i_ang) = bs(n)*sin(ang(i_ang))

            div(i_ang)=y(i_ang)/x(i_ang)

            x1(i_ang) = xs(n) + sqrt(x(i_ang)**2 + y(i_ang)**2) * 
     +           cos(atan(div(i_ang)) - rot(n))
            x2(i_ang) = xs(n) + sqrt(x(i_ang)**2 + y(i_ang)**2) * 
     +           cos(atan(div(i_ang)) - rot(n) - PI)
            y1(i_ang) = ys(n) + sqrt(x(i_ang)**2 + y(i_ang)**2) * 
     +           sin(atan(div(i_ang)) - rot(n))
            y2(i_ang) = ys(n) + sqrt(x(i_ang)**2 + y(i_ang)**2) * 
     +           sin(atan(div(i_ang)) - rot(n) - PI)
            dimen_x1 = dimen_x1+1
            dimen_x2 = dimen_x2+1
         enddo


c initialisation des vecteurs a et b : ils ont deux fois la dimension de x1 
c et x2

         do i=1,dimen_x1
            a(i)=x1(i)
            b(i)=y1(i)
         enddo

c euh je comprends pas trop pouquoi dimen_x2.... car en fait les valeurs vont
c a la suite de dimen_x1.... demande a Leo

         do i=1,dimen_x2
            a(dimen_x2+i)=x2(i)
            b(dimen_x2+i)=y2(i)
         enddo
         
         a(dimen_x1+dimen_x2+1)=x1(1)
         b(dimen_x1+dimen_x2+1)=y1(1)
         
         nb_tot = dimen_x1+dimen_x2+1
         call sm_connect(a,b,nb_tot)

         totox = xs(n) + 1.2 * as(n)
         totoy = ys(n) + 1.2 * bs(n)
         call sm_relocate(totox,totoy)
         ellip = 1. - bs(n)/as(n)
         call sm_expand(0.4)
         write(header,101) ellip
 101     format(f6.4)
         call sm_putlabel(8,header)

      enddo
  
      if(iplot.eq.2)then
         call sm_hardcopy
         write(*,*) "Plot sent to printer"    
      else
         call sm_device(printer)
         call sm_gflush
         call sm_alpha
      endif

      end