c
c------------------------------------------------------------------
c Begin file /home/sting/foster/libplt/mkproj.f
c------------------------------------------------------------------
c
      subroutine mkproj(proj,viewport,perimlat,utc,plat,plon,rot)
c
c Make an ezmap projection (type of projection in proj)
c
c On input:
c   proj = 'CE'   -> make cylindrical equidistant projection
c   proj = 'ST'   -> make stereographic polar projection
c   proj = 'SV'   -> make sat view projection
c   viewport(4) = fractional coords left,right,bottom,top of
c                 desired location of projection
c   perimlat = latitude of perimeter for stereographic proj
c   utc = current ut (used to determine rotation for stereographic)
c
      character*(*) proj
      dimension viewport(4)
      data dgrid/30./
      data ncalls/0/
      save ncalls
c
      ncalls = ncalls+1
      if (ncalls.eq.1) then
        call mapstc('OU','CO')
c       call mapsti('DO',1)
        call mapsti('G1',0)
        call mapsti('G2',2)
        call mapsti('VS',0)
      endif
c
c Cylindrical equidistant:
c
      if (proj.eq.'CE') then
        call mapsti('EL',0)
        call mappos(viewport(1),viewport(2),viewport(3),viewport(4))
        call mapset('MA',0.,0.,0.,0.)
        call mapint
        call maproj(proj,plat,plon,rot)
        call mapint
c
c Stereographic (polar):
c
      elseif (proj.eq.'ST') then
        if (perimlat.lt.-90..or.perimlat.gt.90.) then
          write(6,"('>>> mkproj: bad perimlat=',e12.4)") perimlat
          return
        endif
        plon = 0.
        if (perimlat.le.0.) then
          rrot = utc*15.-180.
          pplat = -90.
        else
          rrot = -utc*15.
          pplat = 90.
        endif
        call mapsti('EL',1)
        call mappos(viewport(1),viewport(2),viewport(3),viewport(4))
        call maproj(proj,pplat,plon,rrot)
        plat = 90.-abs(perimlat)
        call mapset('AN',pplat,pplat,pplat,pplat)
        call mapint
c
c Satellite view projection:
c
      elseif (proj.eq.'SV') then
        call mapsti('EL',1)
        call mapstc('OU','CO')
c       call mapsti('DO',1)
        call mappos(viewport(1),viewport(2),viewport(3),viewport(4))
        call mapstr('SA',rot)
        call mapset('MA',0.,0.,0.,0.)
        call mapstr('GR',dgrid)
        call mapsti('LA',0) ! LA=0 suppresses labeling of meridians and poles
        call mapstr('GD',0.1) ! distance (deg) between pts of grid lines
        call mapint
        call maproj(proj,plat,plon,0.)
        call mapint
c
c Bad proj string:
c
      else
        write(6,"('>>> mkproj: unsupported projection: ',a)") proj
      endif
      return
      end