program triang
      implicit double precision (a-h,o-z)
      integer maxtow,ntower,tottow
c     maxtow is the total number of towers used in the study
c     ntower is the maximum number of towers used for any one location
      parameter (maxtow=100)
      common /param/pi
      double precision pi
      integer itower(maxtow),itowtp(maxtow)
c     tutm are the UTM coordinates of the towers
      double precision tutm(2,maxtow),az(maxtow),coor(2),kappa,vc(2,2),
     1 elarea,area,aztmp(maxtow),sd
      character*40 xeqlin(20)
      integer nargs
      logical plotfg,mle,andrew,huber
      character*40 filein, filout
      character*80 string
      character*80 formt
      data filein/'azimuth.dat'/,filout/'fixes.dat'/,plotfg/.false./,
     1   mle/.false./,huber/.false./,andrew/.true./
c
      pi=datan(1.d0)*(4.d0)
      ntower=3
      formt='(28x,3(i2,f4.0))'
c     set the standard deviation for ALL towers
c     A major problem with the Lenth estimators is that the
c     standard deviation is assumed the same for all towers.
      sd=1.020
c
c   get the execution line commands
c
      call args(nargs,xeqlin)
      do 5 i=1,nargs
         call cvtcas(xeqlin(i))
         if (xeqlin(i)(1:2).eq.'I=') then
            filein=xeqlin(i)(3:)
         else if (xeqlin(i)(1:2).eq.'L=') then
            filout=xeqlin(i)(3:)
         else if (xeqlin(i)(1:3).eq.'SD=') then
            sd=atof(xeqlin(i)(4:))
         else if (xeqlin(i)(1:7).eq.'NTOWER=') then
            ntower=atof(xeqlin(i)(8:))
         else if (xeqlin(i)(1:7).eq.'FORMAT=') then
            formt=xeqlin(i)(8:)
         else if (xeqlin(i).eq.'PLOT') then
            plotfg=.true.
         else if (xeqlin(i).eq.'MLE') then
            mle=.true.
            huber=.false.
            andrew=.false.
         else if (xeqlin(i).eq.'HUBER') then
            mle=.false.
            huber=.true.
            andrew=.false.
         else if (xeqlin(i).eq.'ANDREW') then
            mle=.false.
            huber=.false.
            andrew=.true.
      endif
 5    continue
c
c   read in tower coordinates from the file 'towers.dat'
c
      in=4
      open(unit=in,file='towers.dat',status='old')
      tottow=1
 3    read(in,*,end=4) (tutm(i,tottow),i=1,2)
      tottow=tottow+1
      if (tottow.gt.maxtow) then
         write(*,*) 'Number of tower coordinates exceeds dimensions.'
         stop 'TRIANG bombed'
      endif
      go to 3
 4    tottow=tottow-1
      if (tottow.lt.2) then
         write(*,*) 'Number of tower coordinates is < 2.'
         stop 'TRIANG bombed'
      elseif (tottow.lt.ntower) then
         write(*,*) 'Number of tower coordinates read from file'
         write(*,*) 'TOWERS.DAT is less than number specified'
         write(*,*) 'on execution line.'
         stop 'TRIANG bombed'
      endif
      close(unit=4,status='keep')
c
c   open up file containing azimuths
c
      open(unit=in,file=filein,status='old')
      iout=11
      open(unit=iout,file=filout,status='unknown')
c
c   read input data and process
c
 10   read(in,fmt='(a)',end=30) string
      read(string,formt) (itowtp(i),aztmp(i),i=1,ntower)
      ntwuse=0
      do 15 i=1,ntower
         if (aztmp(i).eq.0.d0) go to 15
         if (itowtp(i).le.0 .or. itowtp(i).gt.maxtow) go to 15
         ntwuse=ntwuse+1
         itower(ntwuse)=itowtp(i)
         az(itower(i))=(90.-abs(aztmp(i)))*(pi/180.)
 15      continue
c     Set ijob to indicate the standard deviation is known.
      ijob=1
      if (ntwuse.ge.2) then
         if (mle) then
         call trimle(tutm,az,itower,ntwuse,coor,sd,kappa,vc,
     1      ijob,ierr)
         else if (huber) then
         call trihub(tutm,az,itower,ntwuse,coor,sd,kappa,vc,
     1      ijob,ierr)
         else
         call triand(tutm,az,itower,ntwuse,coor,sd,kappa,vc,
     1      ijob,ierr)
         endif
c         get confidence ellipse area in hectares
         area=elarea(vc)
c
c         call azplot to plot tower locations and the estimated location
c         on the screen to verify the quality
c
         if (plotfg)
     1       call azplot(string,tutm,itower,az,ntwuse,coor,area)
      else
         coor(1)=0.
         coor(2)=0.
         vc(1,1)=0.
         vc(1,2)=0.
         vc(2,1)=0.
         vc(2,2)=0.
         area=0.
      endif
      write(iout,'(a,2i8.7,f8.5,3f8.1,2i2)')
     1 string(1:lstblk(string)),int(coor(1)+0.5),int(coor(2)+0.5),
     2 area,vc(1,1),vc(1,2),vc(2,2),ntwuse,ierr
      go to 10
 30   continue
      stop 'triang completed normally'
      end