c PROFSPEC - see prof_programs.hlp
c
c  Program plots amplitude, Im part, or Re part vs frequency on semi-log
c  or linear plot from a PROFFT-formatted coefficient file.
c  Tien Grauch, additions by Lin Cordell.
c
      character*50 infile
      character*60 ititle
      character*8 jfmt,jfmt1
      character answer*1,ifmt*56
      dimension isl(5),dxp(2),dyp(2),xp(4),yp(4),ie(10)
      integer kfmt(2),kfmt1(2),jtitle(15)
      equivalence (jfmt,kfmt(1)),(jfmt1,kfmt1(1)),(ititle,jtitle(1))
      data jfmt/'(f8.2)'/,jfmt1/'(g14.6)'/
      icode=0
      isl(1)=2
      iplotr=9
      print803
803   format(' Enter input coefficient file: '$)
      read 804,infile
804   format(a50)
1     print 809
809   format(' 1 - plot amplitude spectrum'/' 2 - plot frequency vs.
     1 Re part', /,' 3 - plot frequency vs. Im part'/,' 4 - plot frequen
     2cy vs. phase angle'/)
      print 808
808   format(' Enter option: ',$)
      read*,icode
      if(icode.ge.1.and.icode.le.4) go to 5
      go to 1
5     open(10,file=infile,status='old',form='formatted',err=298)
      xp(3)=1.5e0
      yp(3)=2.0e0
      xp(2)=0.
      read(10,800) nz,dx
800   format(14x,i1,16x,e16.5)
      read(10,801)
      read(10,801)
801   format(1x)
	print *,'Want program to find min & max for vertical scale? '
	read 888,answer
888	format(a1)
	if(answer.eq.'y'.or.answer.eq.'Y') goto 555
	print*,'enter min & max for vertical scale'
	read*,dyp(1),dyp(2)
	goto 666
555   amax=-1.e38
      amin=1.e38
      frmax=-1.e8
      frmin=1.e38
      fimax=-1.e38
      fimin=1.e38
c  find max and mins of coefs
      i=1
111   read(10,807,end=112) fr,fi,amp
807   format(20x,3e16.5)
      i=i+1
      if(fr.eq.0.0.or.fi.eq.0.0.or.amp.eq.0.0) go to 111
      if(fr.gt.frmax) frmax=fr
      if(fr.lt.frmin) frmin=fr
      if(fi.gt.fimax) fimax=fi
      if(fi.lt.fimin) fimin=fi
      if(amp.gt.amax) amax=amp
      if(amp.lt.amin) amin=amp
      pmax=180
      pmin=-180
      go to 111
112   iend=i
	dyp(2)=0.
666   rewind 10
      read(10,801)
      read(10,801)
      read(10,801)
      dxp(2)=1.1*.5/dx
      if(nz.eq.2) dxp(2)=1.1/dx
      dxp(1)=0.
      call fixscale(dxp(2),dxp(1),xdel)
      go to (700,701,702,703), icode
700   ifmt='(4x,e16.5,32x,e16.5)'
	if(dyp(2).ne.0.) goto 10
      dyp(2)=amax*1.5
      dyp(1)=amin*.8
      go to 10
701   ifmt='(4x,2e16.5)'
	if(dyp(2).ne.0.) goto 10
      dyp(2)=frmax*1.5
      dyp(1)=frmin*.8
      go to 10
702   ifmt='(4x,e16.5,16x,e16.5)'
	if(dyp(2).ne.0.) goto 10
      dyp(2)=fimax*1.5
      dyp(1)=fimin*.8
703   ifmt='(4x,3e16.5)'
	if(dyp(2).ne.0.) goto 10
      dyp(2)=pmax
      dyp(1)=pmin
      go to 10
10    call fixscale(dyp(2),dyp(1),ydel)
	print 815, dxp(1),dxp(2),dyp(1),dyp(2)
815	format(' x values:  minimum = ',g10.4,' maximum = ',g10.4,
     1 /' y values:  minimum = ',g10.4,' maximum = ',g10.4)
      print*,'Enter title (max 60 chars)'
      read 805, ititle
805   format(a)
	if (dyp(1).gt.0.0) then
		print*,'Semi-log or linear plot [s or l]?'
		read 806, answer
		yp(2)=1.0
		if(answer.eq.'l'.or.answer.eq.'L') yp(2)=0.0
	else
		yp(2)=0.0
	endif
15    print 811
811   format(' Plot on screen(cga,ega,vga) or hp plotter(hp): '$)
      read 806, answer
806   format(a1)
      if(answer.ne.'h'.and.answer.ne.'H') go to 25
      iplotr=5
      call pltset(iplotr,xboard,yboard,isl)
      print 810,xboard-1.5,yboard-2.
810   format(' Max size available (in cm) is ',g10.4,' X ',g10.4)
      print*,' Enter vertical axis length in cm'
      read*,yp(1)
20    print*,' Enter horiz axis length in cm'
      read*,xp(1)
      go to 30
25    if(answer.eq.'c'.or.answer.eq.'C') then
      iplotr=8
      else if(answer.eq.'e'.or.answer.eq.'E') then
      iplotr=9
      else if(answer.eq.'v'.or.answer.eq.'V') then
      iplotr=10
      else
      go to 15
      endif
      call pltset(iplotr,xboard,yboard,isl)
      xp(1)=xboard-xp(3)-1
      yp(1)=yboard-yp(3)-1
30    nopts=4
      xp(4)=xp(1)+xp(3)+1
      yp(4)=yp(1)+yp(3)+1
      call scale(dxp,dyp,xp,yp,nopts,ier)
40    continue
      if(icode.ne.4)then
        read(10,ifmt,end=50) x,y
      else
        read(10,ifmt,end=50)x,fr,fi
        if(fr.eq.0.and.fi.eq.0)then
          print*,'Indeterminate phase at',x
          phase=0.
        else
          phase=atan2(-fi,fr)
          phase=(phase*180.)/3.1415927
        endif
        y=phase
      endif
c  if y <= 0 cannot plot on log scale, so plot off the scale
      if(yp(2).ne.0.0.and.y.le.0.0) y=1.e+30
      call vchar(x,y,4,1,1,.16,0.,0.,0.)
      go to 40
50    call xaxis(dxp,dyp,xp,xdel,5,.2,kfmt,8)
	if (yp(2).ne.0.0) then
	      call yaxsl(dyp,dxp,yp,.2,kfmt1,14)
	else
	      call yaxis(dyp,dxp,yp,ydel,5,.2,kfmt,8)
	endif
60    xd=3.5
      call vchar(xd,.5,jtitle,60,3,.25,0.,0.,0.)
      call neatl
      ie(2)=2
      call endpt(ie)
      go to 300
298   print*,'Input file not found or not binary.'
300   close(10)
       stop
       end
      subroutine fixscale(xmax,xmin,del)
c
c  For plot_sys programs to find parameter del for call to xaxis or
c  yaxis knowing end points of axis
c
      power=-1.0
      constant=.5
      idel=alog10(xmax-xmin)
      tmp=float(1-idel)
      ten=10.e0**tmp
      quantity=ten*(xmax-xmin)
      if(quantity.gt.30.0) constant=1.0
      if(quantity.lt.15.0) constant=.2
      del=ten**power
      del=del*constant
      return
      end
