      real*8 PI,RAD
      parameter (PI=3.141592653589793, RAD=180.0/PI)
      real lmd,pshi,x,y

      character name(1000)*10,cmnt(1000)*13,fname1*90, fname2*90
      real alpha(1000),delta(1000),gl(1000),gb(1000),intens(1000)
      character strnam*10,nname*15,dummy*80
      character sname*60


      parameter (nlim=1800)
      real*4 data(nlim,3)
      character aname(250)*15
      real adata(250,nlim,3),ax(250),ay(250)

      fname1=
     &   "/home/kitamoto/Asm/Asm/analysis_tools/asmmap/asmsource"
      fname2=
     &   "/home/kitamoto/Asm/Asm/analysis_tools/asmmap/source"

      call pgbegin(0,'?',1,1)
      call pgenv(-320.,320.,-160.,160.,0,-2)

      call waku

 2000 continue
      open(7, file=fname1, status='old')
      print *,"opened",fname
      read(7,'(//)')
      idata=1
 2100 continue
      print *,"idata=",idata
      read(7,'(i3,1x,a10,1x,2f8.4,1x,f9.4,f8.4,2x,f7.2,a13,8x)',
     &     end=2200)
     &     id,name(idata),alpha(idata),delta(idata),
     &     gl(idata),gb(idata),intens(idata),cmnt(idata)
      idata=idata+1
      goto 2100
 2200 continue
      close(7)

      ndata=0
      call pgsch(0.8)

*      open(8,file='/home/kitamoto/Asm/Asm/analysis_tools/asmmap/source',
      open(8,file=fname2,
     &     status='old')
 3000 CONTINUE
      ndata=ndata+1
      READ(8,*,END=3200) STRNAM,nname
      read(8,*,end=3200) dummy
      aname(ndata)=nname
       print *,"dummy=", dummy(1:20)
       print *," fname2=", fname2

c      read(dummy,'(12x15a)') sname
*      write(sname,'(''~egoshi/ASM.database/'',a)') dummy 
*      write(sname,'(''/export/home/kitamoto/Asm/asmmap/ASM.database/'',a)') dummy
*       print *,"dummy  =  ", dummy
*       print *,"STARNAM , nname =  ", STANAM, nname
      write(sname,'(''/home/kitamoto/Asm/Asm/data/'',a)') dummy(1:20)
      call readd2(sname,data,iret)
      write(*,*) dummy(13:27),iret
 3301 do 3300 i33=1,nlim
 3351    do 3350 i335=1,3
            adata(ndata,i33,i335)=data(i33,i335)
 3350    continue
 3300 continue
 3101 do 3100 i31=1,idata
         if (strnam.eq.name(i31)) then
            lmd=gl(i31)/RAD
            pshi=gb(i31)/RAD
            call etv(lmd,pshi,x,y)
            x=-x
            ax(ndata)=x
            ay(ndata)=y
            rinten=adata(ndata,nlim,3)
            if(rinten.lt.0.01) then
               call pgsci(15)
            elseif(rinten.lt.0.5) then
               call pgsci(1)
            elseif(rinten.lt.1.0) then
               call pgsci(5)
            elseif(rinten.lt.1.5) then 
               call pgsci(8)
            else
               call pgsci(7)
            endif
            if (rinten.lt.0.01) rinten=0.01
            call pgsch(log10(rinten*200))
            call pgpoint(1,x,y,17)
         endif
 3100 continue
      goto 3000
 3200 continue

      call pgask(.false.)
 5001 do 5000 i5=70,nlim
         call pgpage
         call pgsci(1)
         call waku
         write(*,*) i5
 5501    do 5500 i55=1,ndata
            rinten=adata(i55,i5,3)
            if(rinten.lt.0.01) then
               call pgsci(15)
            elseif(rinten.lt.0.5) then
               call pgsci(1)
            elseif(rinten.lt.1.0) then
               call pgsci(5)
            elseif(rinten.lt.1.5) then 
               call pgsci(8)
            else
               call pgsci(7)
            endif
            if (rinten.lt.0.01) rinten=0.01
            call pgsch(log10(rinten*200))
            call pgpoint(1,ax(i55),ay(i55),17)
 5500    continue
 5901    do 5900 i59=1,500000000
 5900    continue
 5000 continue

      call pgsci(1)
      call pgsch(1.5)
      call pgptext(0.0,-180.0,0.0,0.5,'Galactic Longitude')
      call pgptext(-340.0,0.0,90.0,0.5,'Galactic Latitude')
      call pgsch(2.0)
      call pgptext(0.0,180.0,0.0,0.5,'\\frGINGA ASM MAP')
      read(*,*)
      call pgend
      stop
      end

c====================================================================
      subroutine etv(lmd,pshi,x,y)
c--------------------------------------------------------------------
      implicit double precision  ( a - h, o - z )
      real lmd,pshi,x,y
      real*8 PI, ETVRAD
      parameter (PI=3.141592653589793d0, ETVRAD=100.0d0)
      cosdlt=cos(pshi)*cos(0.5d0*lmd)
c      sindlt=sqrt(1.00001-cosdlt*cosdlt)
      sindlt=sqrt(1.0000001d0-cosdlt*cosdlt)
c      if (abs(cosdlt) .lt. 0.00001) then
      if (abs(cosdlt) .lt. 1.0d-7) then
         dlt=PI/2.0d0
         goto 1000
      endif
      dlt=atan(sindlt/cosdlt)
 1000 continue
      cosv=sin(pshi)/sin(dlt)
c      sinv=sqrt(1.00001-cosv*cosv)
      sinv=sqrt(1.0000001d0-cosv*cosv)
      x=2.0d0*ETVRAD*dlt*sinv* sign(1.0,lmd)
      y=ETVRAD*dlt*cosv
      return
      end


c====================================================================
      subroutine waku
c--------------------------------------------------------------------
      real*8 PI,RAD
      integer ph, ph2
      parameter (PI=3.141592653589793, RAD=180.0/PI)
      real lmd,pshi,x,y
      real xr(2500),yr(2500)

      ndata=0
      inv=-1
 1001 do 1000 lm=-180,180,15
         inv=inv*(-1)
         lmd=lm/RAD
 1101    do 1100 ph=-90,90,5
            pshi=ph/RAD
            pshi=pshi*inv
            call etv(lmd,pshi,x,y)
            ndata=ndata+1
            xr(ndata)=x
            yr(ndata)=y
 1100    continue
 1000 continue

 1401 do 1400 ph=80,-80,-10
         inv=inv*(-1)
         pshi=ph/RAD
 1551    do 1550 ph2=ph+10,ph,-5
            pshi=ph2/RAD
            call etv(lmd,pshi,x,y)
            ndata=ndata+1
            xr(ndata)=x
            yr(ndata)=y
 1550    continue
 1501    do 1500 lm=-180,180,5
            lmd=lm/RAD
            lmd=lmd*inv
            call etv(lmd,pshi,x,y)
            ndata=ndata+1
            xr(ndata)=x
            yr(ndata)=y
 1500    continue
 1400 continue
      call pgline(ndata,xr,yr)

      return
      end

c====================================================================
      subroutine readd2(fname,data,iret)
c--------------------------------------------------------------------
      parameter (nlim=1800)
      real*4 data(nlim,3),edata(nlim,3)
      integer nchk(nlim)
      character cmnt*20,fname*60
      character dummy*80
      logical exflag
*------- clear data ------------
 2001 do 2000 i21=1,nlim
 2101    do 2100 i22=1,3
	    data(i21,i22)=0.0
 2100    continue
 2000 continue
      iret=0
*------- inquire file ----------
      inquire(file=fname,exist=exflag)
      if(exflag)then
*------- open file -------------
         open(7,file=fname,status='old')
*------- read comment ----------
         read(7,'(a20)') cmnt
         ndata=0
 1000    continue
*------- read data -------------
         read(7,'(a80)',end=999) dummy
         read(dummy,'(i4)') no
         if(no.eq.0) goto 1000
         read(dummy,'(4x,i4,4f8.4,2f8.3)') ich,d,e1,e2,ed,da,el
         if(ich.eq.108) then
            mch=1
         elseif(ich.eq.916) then
            mch=2
         elseif(ich.eq.116) then
            mch=3
         else
            goto 1000
         endif
         iday=da
         data(iday,mch)=d
         edata(iday,mch)=ed
*------- next data -------------
         goto 1000
 999     continue
         close(7)
ce------ Check bad data --------
 1401    do 1400 i14=1,nlim
            nchk(i14)=1
 1411       do 1410 i141=1,2
               if(data(i14,i141).le.-0.3) then
                  edata(i14,i141)=0
                  goto 1410
               elseif(edata(i14,i141).ge.0.5) then
                  edata(i14,i141)=0
                  goto 1410
               elseif(edata(i14,i141).le.2e-3) then
                  edata(i14,i141)=0
                  goto 1410
               elseif(data(i14,i141).le.0) then
                  zero = data(i14,i141) + edata(i14,i141)*3
                  if(zero.lt.0) then
                     edata(i14,i141)=0
                     goto 1410
                  endif
               endif
 1410       continue
 1400    continue
*------- Interpolate data -------------
 1531    do 1530 i153=2,nlim
            if(edata(i153,1).eq.0 .or. edata(i153,2).eq.0) then
               data(i153,1) = data(i153-1,1)
               data(i153,2) = data(i153-1,2)
            endif
 1530    continue
*------- for all channel -------
 1301    do 1300 i13=1,nlim
            if(data(i13,3).eq.0.0) then
               data(i13,3)=data(i13,1)+data(i13,2)
            endif
 1300    continue
*------- return ----------------
      else
         write(*,*) 'File cannot be found ! Input once again!'
         write(*,*) fname
         iret=9
      endif
 990  continue
      return
      end
