*=================================================================
*
*   Program for analysis of the ASM data
*
*				By Hiroaki Watanabe
*
*				1992.11.20
*
*=================================================================
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	integer	  rch(8)
	integer   ipage
ce	character cmnt*20,fname*15,fname1*25
	character cmnt*20,fname*15,fname1*60
	character sname*20
ce
        logical exflag
*------- Program start !! ---------
	write(*,*) '*============================================*'
	write(*,*) '*				              *'
	write(*,*) '*          Welcome to ASM world !!           *'
	write(*,*) '*				              *'
	write(*,*) '*    This program  is  very  useful for  an  *'
	write(*,*) '*  analysis for the ASM data.  The original  *'
	write(*,*) '*  data  are  stored   in  the  side  B  of  *'
	write(*,*) '*  Magnetic  Optical  Disk.  Please confirm  *'
	write(*,*) '*  that the data for each source  which you  *'
	write(*,*) '*  would like to analyze  have been made by  *'
	write(*,*) '*  another program. If some troubles happen, *'
	write(*,*) '*  would you make contact with the author ?  *'
	write(*,*) '*				              *'
	write(*,*) '*        by H.Watanabe on November 20, 1992. *'
	write(*,*) '*				              *'
	write(*,*) '*============================================*'
	write(*,*)
*------- Set PGPLOT page ----------
	ipage=1
*------- Input source name --------
 1000	continue
	write(*,*) 'Please input source file name-- !'
	read(*,'(a)') fname
	write(*,*) fname
	if(fname.eq.'q') goto 999
ce---  	write(fname1,'(''data/'',a)') fname
*    	write(fname1,'(''/export/home/kitamoto/Asm/data/'',a)') fname 
        write(fname1,'(''/home/kitamoto/Asm/Asm/data/'',a)') fname

	print *,"fname1=",fname1

        inquire(file=fname1,exist=exflag)
        if(.not.exflag)then
           write(fname1,'(''/home/kitamoto/Asm/data/'',a)')
     &          fname
           print *,"fname1=",fname1
        endif
*------- Read data ----------------
 	call readd(fname1,day,data,edata,ndata,cmnt,iret)

*	     print *,"just after readd ndata=", ndata
*	     print *,(j, day(j), data(j, 17), edata(j,17),j=1,10)

	if(iret.eq.9) goto 1000
	write(*,*) 'Input Nick Name ?'
	read(*,'(a)') sname
	if(sname.eq.'q') then
	  cmnt(1:20)=cmnt(2:11)//'          '
	else
	  cmnt(1:20)=sname
	endif
	write(*,*)
*------- Menu ---------------------
 2000	continue
	write(*,*) '*==============      Menu      ==============*'
	write(*,*)
	write(*,*) '	 1. Light curve (1-6 keV & 6-20 keV)'
	write(*,*) '	 2. Light curve (1-20 keV)'
	write(*,*) '	 3. Light curve (each channel)'
	write(*,*) '	 4. Hardness (Sofeness) Ratio'
	write(*,*) '	 5. All 16-channels plot'
	write(*,*) '	 6. Color - Color Diagram'
	write(*,*) '	 7. Intensity - Color Diagram'
	write(*,*) '	 8. Change the source file name'
	write(*,*) '	 9. Intensity Distridution'
	write(*,*) '	10. Spectrum'
	write(*,*) '	11. Mac file'
	write(*,*) '	99. End'
	write(*,*)
	write(*,199) cmnt(1:20)
  199	format(12x,'Source name = ',a20)
	write(*,*)
	write(*,*) '	    Please select the number!'
*------- Input number -------------
	read(*,'(i2)') num
	write(*,*)
*------- Number 1 -----------------
	if(num.eq.1) then
	  xmin=0.0
	  xmax=0.0
	  ymin=0.0
	  ymax=0.0
	  write(*,*) 'Linear (1) or Log (2) ?'
	  read(*,'(i2)') ill
	  write(*,*)
 	  if(ill.eq.1) then
 	  call dplot3
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  else
 	  call dplot4
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  endif
*------- Change Scale -------------
  130	  continue
	  write(*,*)
	  write(*,*) 'Change the scales (y-1/n-0) ?'
	  read(*,'(i2)') iyn
	  write(*,*)
	  if(iyn.ne.1) goto 120
	  call scalet(xmin,xmax,ymin,ymax)
	  if(ill.eq.1) then
 	  call dplot3
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  else
 	  call dplot4
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  endif
	  goto 130
  120	continue
	endif
*------- Number 2 -----------------
	if(num.eq.2) then
	  nch=19
	  xmin=0.0
	  xmax=0.0
	  ymin=0.0
	  ymax=0.0
	  write(*,*) 'Linear (1) or Log (2) ?'
	  read(*,'(i2)') ill
	  write(*,*)
	  if(ill.eq.1) then
	  call dplot1
     &    (nch,day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  else
 	  call dplot2
     &    (nch,day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  endif
*------- Change Scale -------------
  230	  continue
	  write(*,*)
	  write(*,*) 'Change the scales (y-1/n-0) ?'
	  read(*,'(i2)') iyn
	  if(iyn.ne.1) goto 220
	  call scalet(xmin,xmax,ymin,ymax)
	  if(ill.eq.1) then
	  call dplot1
     &    (nch,day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  else
 	  call dplot2
     &    (nch,day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  endif
	  goto 230
  220	  continue
	endif
*------- Number 3 -----------------
	if(num.eq.3) then
	  call lcurve(day,data,edata,ndata,cmnt,ipage)
	endif
*------- Number 4 -----------------
	if(num.eq.4) then
	  call hratio(day,data,edata,ndata,cmnt,ipage)
	endif
*------- Number 5 -----------------
	if(num.eq.5) then
	  xmin=0.0
	  xmax=0.0
	  ymin=0.0
	  ymax=0.0
	  write(*,*) 'Linear (1) or Log (2) ?'
	  read(*,'(i2)') ill
	  write(*,*)
 	  if(ill.eq.1) then
 	  call dplot10
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  else
 	  call dplot11
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  endif
*------- Change Scale -------------
 1130	  continue
	  write(*,*)
	  write(*,*) 'Change the scales (y-1/n-0) ?'
	  read(*,'(i2)') iyn
	  write(*,*)
	  if(iyn.ne.1) goto 1120
	  call scalet(xmin,xmax,ymin,ymax)
	  if(ill.eq.1) then
 	  call dplot10
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  else
 	  call dplot11
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,ymin,ymax,ipage)
	  endif
	  goto 1130
 1120	continue
	endif
*------- Number 6 -----------------
	if(num.eq.6) then
	  write(*,*) '*--------- Color - Color Diagram --------*'
	  write(*,*)
*------- Input channel ---------
  640	  continue
*------- Yoko jiku -------------
	  write(*,*) 'Input Color (A/B) channel (A1-A2/B1-B2)'
	  write(*,*)
	  write(*,*) 'A1 ?'
	  read(*,'(i2)') rch(1)
	  write(*,*) 'A2 ?'
	  read(*,'(i2)') rch(2)
	  write(*,*) 'B1 ?'
	  read(*,'(i2)') rch(3)
	  write(*,*) 'B2 ?'
	  read(*,'(i2)') rch(4)
	  write(*,*)
*------- Judgement -------------
	  if(rch(1).gt.rch(2).or.rch(3).gt.rch(4)) then
	    write(*,*) 'Input once again !!'
	    write(*,*)
	    goto 640
	  endif
	  if(rch(1).gt.16.or.rch(2).gt.16) goto 640
	  if(rch(3).gt.16.or.rch(4).gt.16) goto 640
	  write(*,611) (rch(i1),i1=1,4)
  611	  format(' Color (A/B) = ',i3,'-',i3,'/',i3,'-',i3)
	  write(*,*)
*------- Tate jiku -------------
	  write(*,*) 'Input Color (C/D) channel (C1-C2/D1-D2)'
	  write(*,*)
	  write(*,*) 'C1 ?'
	  read(*,'(i2)') rch(5)
	  write(*,*) 'C2 ?'
	  read(*,'(i2)') rch(6)
	  write(*,*) 'D1 ?'
	  read(*,'(i2)') rch(7)
	  write(*,*) 'D2 ?'
	  read(*,'(i2)') rch(8)
	  write(*,*)
*------- Judgement -------------
	  if(rch(5).gt.rch(6).or.rch(7).gt.rch(8)) then
	    write(*,*) 'Input once again !!'
	    write(*,*)
	    goto 640
	  endif
	  if(rch(5).gt.16.or.rch(6).gt.16) goto 640
	  if(rch(7).gt.16.or.rch(8).gt.16) goto 640
	  write(*,612) (rch(i2),i2=5,8)
  612	  format(' Color (C/D) = ',i3,'-',i3,'/',i3,'-',i3)
	  write(*,*)
*------- Plot ------------------
 	  call color(rch,day,data,edata,ndata,cmnt,ipage)
*------- Next color ------------
          write(*,*) 'Another one (y-1/n-0) ?'
	  read(*,'(i2)') iyn
	  write(*,*)
	  if(iyn.ne.1) goto 650
	  goto 640
  650	  continue
	endif
*------- Number 7 -----------------
	if(num.eq.7) then
*------- Input channel ---------
	  write(*,*) '*------- Intensity - Color Diagram -------*'
	  write(*,*)
  940	  continue
*------- Input A ---------------
  941	  continue
	  write(*,*) 'Input Intensity channel (A1-A2)'
	  write(*,*)
	  write(*,*) 'A1 ?'
	  read(*,'(i2)') rch(1)
	  write(*,*) 'A2 ?'
	  read(*,'(i2)') rch(2)
	  write(*,*)
*------- Judgement -------------
	  if(rch(1).gt.16.or.rch(2).gt.16) goto 941
	  if(rch(1).gt.rch(2)) goto 941
	  write(*,944) (rch(i4),i4=1,2)
  944	  format(' Intensity channel ',i3,'-',i3)
	  write(*,*)
*------- Tate jiku -------------
  942	  continue
	  write(*,*) 'Input Color channel (B1-B2/C1-C2)'
	  write(*,*)
	  write(*,*) 'B1 ?'
	  read(*,'(i2)') rch(3)
	  write(*,*) 'B2 ?'
	  read(*,'(i2)') rch(4)
	  write(*,*)
	  write(*,*) 'C1 ?'
	  read(*,'(i2)') rch(5)
	  write(*,*) 'C2 ?'
	  read(*,'(i2)') rch(6)
	  write(*,*)
*------- Judgement -------------
	  if(rch(3).gt.16.or.rch(4).gt.16) goto 942
	  if(rch(5).gt.16.or.rch(6).gt.16) goto 942
	  if(rch(3).gt.rch(4)) goto 942
	  if(rch(5).gt.rch(6)) goto 942
	  write(*,943) (rch(i3),i3=3,6)
  943	  format(' Color channel ',i3,'-',i3,' /',i3,'-',i3)
	  write(*,*)
*------- Plot ------------------
 	  call icolor(rch,day,data,edata,ndata,cmnt,ipage)
*------- Next color ------------
          write(*,*) 'Another one (y-1/n-0) ?'
	  read(*,'(i2)') iyn
	  write(*,*)
	  if(iyn.ne.1) goto 950
	  goto 940
  950	  continue
	endif
*------- Number 8 -----------------
	if(num.eq.8) goto 1000 
*------- Number 9 -----------------
	if(num.eq.9) then
	  call dplot8
     &    (day,data,edata,ndata
     &    ,cmnt,xmin,xmax,nbin,ipage)
	endif
*------- Number 10 ----------------
	if(num.eq.10) then
	  call spec(day,data,edata,ndata,cmnt,ipage)
	endif
*------- Number 11 ----------------
	if(num.eq.11) then
	  if(fname1.eq.'data/'.or.ndata.le.0) then
	    write(*,*) 'No data in the buffer.'
	    write(*,*)
	    goto 7000
	  else
	     print *,"ndata=", ndata
*	     print *,(j, day(j), data(j, 17), edata(j,17),j=1,10)
	    call wrtfil(fname,day,data,edata,ndata,cmnt)
	  endif
	endif	  
*------- Number 99 ----------------
	if(num.eq.99) goto 999
*------- Goto menu ----------------
 7000	continue
	goto 2000
*------- Program end --------------
  999	continue
	if(ipage.ne.1) then
	  call pgend
	  write(*,*)
	  write(*,*)
	endif
	write(*,*) ' ****** See you again. Good-bye! ******'
	write(*,*)
	stop
	end	
******************************************************************
	subroutine dplot1
     & (nch,day,data,edata,ndata,cmnt,xmin,xmax,ymin,ymax,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20
	character text*14
	integer   ipage
*------- plot channel ----------
	write(text,'(''1-20 keV  '')')
*------- scale -----------------
	if(xmax.eq.0.0) xmax=1800.0
	if(ymax.eq.0.0) then
	  ymax1=data(1,nch)
 2001	  do 2000 i2=2,ndata
	    if(ymax1.lt.data(i2,nch)) ymax1=data(i2,nch)
 2000	  continue
	  if(ymax1.le.0.0) ymax1=1.0
	else
	  ymax1=ymax
	endif
*------ Light Curve ------------
	if(ipage.eq.1) call pgbegin(0,'?',1,1)
	call pgsch(0.9)
*------ Waku wo kaku -----------
	if(ymax.le.0.0) then
	if(ymax1.le.1.0) then
	  ymax1=ymax1*100000.
	  iflag=1
	endif
	call sutble(ymax1,ymax)
	if(iflag.eq.1) ymax=ymax/100000.
	endif
	call pgsci(4)
	call pgscf(1)
	call pgenv(xmin,xmax,ymin,ymax,0,0)
*------ plot comment -----------
	call pgsch(0.97)
	call pgsci(5)
	call pgscf(2)
ce	call pglabel('Day of 1987','Flux (counts/sec/cm**2)',
	call pglabel('Day of 1987','Flux (counts/sec/cm\\u2\\d)',
     &	             'ASM Light Curve of '//cmnt(1:20))
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.85)
 2201	do 2200 i21=1,ndata
	  if(edata(i21,nch).le.0.0) goto 2200
	  x0=day(i21)
	  y0=data(i21,nch)
	  y1=data(i21,nch)-edata(i21,nch)
	  y2=data(i21,nch)+edata(i21,nch)
	  call pgsci(7)
	  call pgpoint(1,x0,y0,15)
	  call pgerry(1,x0,y1,y2,0.0)
 2200	continue
*------ title ------------------
	call pgsci(8)
	call pgscf(3)
	call pgsch(0.9)
	x1=xmin+(xmax-xmin)*0.73
	y1=ymin+(ymax-ymin)*0.87
	call pgtext(x1,y1,text//char(15))
	call pgsch(0.8)
	call pgiden
*------ end --------------------
  970	continue
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
************************************************************
	subroutine dplot2
     & (nch,day,data,edata,ndata,cmnt,xmin,xmax,ymin,ymax,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20
	character text*14
	integer   ipage
*------- plot channel ----------
	write(text,'(''1-20 keV  '')')
*------- scale -----------------
	if(xmax.eq.0.0) xmax=1800.0
	if(ymax.le.0.0) then
	  ymax1=data(1,nch)
 2001	  do 2000 i2=2,ndata
	    if(ymax1.lt.data(i2,nch)) ymax1=data(i2,nch)
 2000	  continue
	  if(ymax1.le.0.0) ymax1=1.0
	else
	  ymax1=ymax
	endif
*------ Light Curve ------------
	if(ipage.eq.1) call pgbegin(0,'?',1,1)
	call pgsch(0.9)
*------ Kirino ii suuji --------
	iyord=3
	if(ymax.le.0.0) then
	  iloy=int(alog10(ymax1))+1
	  ylog1=float(iloy-1)-alog10(ymax1)
	  if(ylog1.gt.0.0) iloy=iloy-1
	  ymax=iloy
	else
	  ymax=alog10(ymax1)
	  iloy=int(ymax)+1
	endif
	if(ymin.le.0.0) then
	  ymin0=10.**(iloy-iyord)
	  ymin=iloy-iyord
	else
	  ymin0=ymin
	  ymin=alog10(ymin)
	endif
*------ Waku wo kaku -----------
	call pgsci(4)
	call pgscf(1)
	call pgenv(xmin,xmax,ymin,ymax,0,20)
*------ plot comment -----------
	call pgsch(0.97)
	call pgsci(5)
	call pgscf(2)
ce	call pglabel('Day of 1987','Flux (counts/sec/cm**2)',
	call pglabel('Day of 1987','Flux (counts/sec/cm\\u2\\d)',
     &	             'ASM Light Curve of '//cmnt(1:20))
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.85)
 2201	do 2200 i21=1,ndata
	  x0=day(i21)
	  if(data(i21,nch).le.0.0) then
	    y0=alog10(ymin0*0.1)
	  else
	    y0=alog10(amax1(data(i21,nch),ymin0))
	  endif
	  y1=alog10(amax1(data(i21,nch)-edata(i21,nch),ymin0))
	  y2=alog10(amax1(data(i21,nch)+edata(i21,nch),ymin0))
	  call pgsci(7)
	  call pgpoint(1,x0,y0,15)
	  call pgerry(1,x0,y1,y2,0.0)
 2200	continue
*------ title ------------------
	call pgsci(8)
	call pgscf(3)
	call pgsch(0.9)
	x1=xmin+(xmax-xmin)*0.73
	y1=ymin+(ymax-ymin)*0.87
	call pgtext(x1,y1,text//char(15))
	call pgsch(0.8)
	call pgiden
*------ end --------------------
  970	continue
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
************************************************************
	subroutine sutble(a,b)
	b=0
	if(a.le.0.0) return
	x=alog10(a)
	i=x
	if(a.lt.1.) i=i-1
	x=x-i
	x=10.**x
	b=10
	if(x.le.8.0) b=8
	if(x.le.5.0) b=5
	if(x.le.4.0) b=4
	if(x.le.2.0) b=2
	b=b*10.**float(i)
	return
	end
************************************************************
	subroutine scalet(xmin,xmax,ymin,ymax)
	real xmin,xmax,ymin,ymax
	write(*,*)
     & 'If you set 0.0, the value will be set automatically.'
	write(*,*) 'X minimum value?'
	read(*,*) xmin
	write(*,*) 'X maximum value?'
	read(*,*) xmax
	write(*,*) 'Y minimum value?'
	read(*,*) ymin
	write(*,*) 'Y maximum value?'
	read(*,*) ymax
	write(*,*)
	return
	end	
************************************************************
	subroutine dplot3
     &   (day,data,edata,ndata,cmnt,xmin,xmax,ymin,ymax,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20
	character text1*10,text2*10
	integer   ipage
*------- plot channel ---------
	write(text1,'('' 1-6 keV  '')')
	write(text2,'(''6-20 keV  '')')
*------- Set scale ------------
	if(xmax.le.0.) xmax=1800.0
	if(ymax.le.0.0) then
	  ymax1=data(1,17)
 2001	  do 2000 i2=2,ndata
 3001	    do 3000 i3=17,18
	      if(ymax1.lt.data(i2,i3)) ymax1=data(i2,i3)
 3000	    continue
 2000	  continue
	  if(ymax1.le.0.0) ymax1=1.0
	else
	  ymax1=ymax
	endif
*------ Light Curve -----------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
	call pgsch(0.9)
*------ kirino ii suuji -------
	if(ymax.le.0.0) then
	if(ymax1.le.1.0) then
	  ymax1=ymax1*100000.
	  iflag=1
	endif
	call sutble(ymax1,ymax)
	if(iflag.eq.1) ymax=ymax/100000.
	endif
*------ Plot -----------------
	isym=1
	call pgvport(0.1,0.9,0.545,0.93)
  700 	continue
*------ Waku wo kaku -----------
	call pgsci(3)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin,ymax)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgsch(0.95)
	call pgsci(8)
	call pgscf(2)
	if(isym.eq.1)then
ce	  call pglabel(' ','Flux (counts/sec/cm**2)',
	  call pglabel(' ','Flux (counts/sec/cm\\u2\\d)',
     &	               'ASM Light Curve of '//cmnt(1:20))
	else
ce	  call pglabel('Day of 1987','Flux (counts/sec/cm**2)',' ')
	  call pglabel('Day of 1987','Flux (counts/sec/cm\\u2\\d)',' ')
	endif
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.75)
 2201	do 2200 i21=1,ndata
	x0=day(i21)
	if(isym.eq.1) then
	  if(edata(i21,17).le.0.0) goto 2200
	  y0=data(i21,17)
	  y1=data(i21,17)-edata(i21,17)
	  y2=data(i21,17)+edata(i21,17)
	  call pgsci(7)
	  call pgpoint(1,x0,y0,15)
	else
	  if(edata(i21,18).le.0.0) goto 2200
	  y0=data(i21,18)
	  y1=data(i21,18)-edata(i21,18)
	  y2=data(i21,18)+edata(i21,18)
	  call pgsci(7)
	  call pgpoint(1,x0,y0,18)
	endif
	call pgerry(1,x0,y1,y2,0.0)
 2200	continue
*------ title ------------------
	call pgsci(5)
	call pgscf(3)
	call pgsch(0.9)
	x1=xmin+(xmax-xmin)*0.76
	y1=ymin+(ymax-ymin)*0.85
	if(isym.eq.1) then
	  call pgtext(x1,y1,text1//char(15))
	else
	  call pgtext(x1,y1,text2//char(18))
	endif
*------ Next energy Channel ----
	isym=isym+1
	if(isym.gt.2) goto 970
	call pgvport(0.1,0.9,0.08,0.475)
	goto 700
*------ end --------------------
  970	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
************************************************************
	subroutine dplot4
     &   (day,data,edata,ndata,cmnt,xmin,xmax,ymin,ymax,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20
	character text1*10,text2*10
	integer   ipage
*------- plot channel ---------
	write(text1,'('' 1-6 keV  '')')
	write(text2,'(''6-20 keV  '')')
*------- Set scale ------------
	if(xmax.le.0.) xmax=1800.0
	if(ymax.le.0.0) then
	  ymax1=data(1,17)
 2001	  do 2000 i2=2,ndata
 3001	    do 3000 i3=17,18
	      if(ymax1.lt.data(i2,i3)) ymax1=data(i2,i3)
 3000	    continue
 2000	  continue
	  if(ymax1.le.0.0) ymax1=1.0
	else
	  ymax1=ymax
	endif
*------ Light Curve -----------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
	call pgsch(0.9)
*------ Kirino ii suuji --------
	iyord=3
	if(ymax.le.0.0) then
	  iloy=int(alog10(ymax1))+1
	  ylog1=float(iloy-1)-alog10(ymax1)
	  if(ylog1.gt.0.0) iloy=iloy-1
	  ymax=iloy
	else
	  ymax=alog10(ymax1)
	  iloy=int(ymax)+1
	endif
	if(ymin.le.0.0) then
	  ymin0=10.**(iloy-iyord)
	  ymin=iloy-iyord
	else
	  ymin0=ymin
	  ymin=alog10(ymin)
	endif
*------ Plot -----------------
	isym=1
	call pgvport(0.1,0.9,0.545,0.93)
  700 	continue
*------ Waku wo kaku -----------
	call pgsci(3)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin,ymax)
	call pgbox('BCTSN',0.0,0,'BCTSNL',0.0,0)
*------ plot comment -----------
	call pgsch(0.95)
	call pgsci(8)
	call pgscf(2)
	if(isym.eq.1)then
ce	  call pglabel(' ','Flux (counts/sec/cm**2)',
	  call pglabel(' ','Flux (counts/sec/cm\\u2\\d)',
     &	               'ASM Light Curve of '//cmnt(1:20))
	else
ce	  call pglabel('Day of 1987','Flux (counts/sec/cm**2)',' ')
	  call pglabel('Day of 1987','Flux (counts/sec/cm\\u2\\d)',' ')
	endif
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.75)
 2201	do 2200 i21=1,ndata
	x0=day(i21)
	if(isym.eq.1) then
	  if(data(i21,17).le.0.0) then
	    y0=alog10(ymin0*0.1)
	  else
	    y0=alog10(amax1(data(i21,17),ymin0))
	  endif
	  y1=alog10(amax1(data(i21,17)-edata(i21,17),ymin0))
	  y2=alog10(amax1(data(i21,17)+edata(i21,17),ymin0))
	  call pgsci(7)
	  call pgpoint(1,x0,y0,15)
	else
	  if(data(i21,18).le.0.0) goto 2200
	  y0=alog10(amax1(data(i21,18),ymin0))
	  y1=alog10(amax1(data(i21,18)-edata(i21,18),ymin0))
	  y2=alog10(amax1(data(i21,18)+edata(i21,18),ymin0))
	  call pgsci(7)
	  call pgpoint(1,x0,y0,18)
	endif
	call pgerry(1,x0,y1,y2,0.0)
 2200	continue
*------ title ------------------
	call pgsci(5)
	call pgscf(3)
	call pgsch(0.9)
	x1=xmin+(xmax-xmin)*0.76
	y1=ymin+(ymax-ymin)*0.85
	if(isym.eq.1) then
	  call pgtext(x1,y1,text1//char(15))
	else
	  call pgtext(x1,y1,text2//char(18))
	endif
*------ Next energy Channel ----
	isym=isym+1
	if(isym.gt.2) goto 970
	call pgvport(0.1,0.9,0.08,0.475)
	goto 700
*------ end --------------------
  970	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
********************************************************
	subroutine hratio(day,data,edata,ndata,cmnt,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
ce---	real*4    engch1(16)/1.000,1.674,2.249,2.860,
*    &		             3.479,4.100,4.750,5.338,
*    &		             5.920,7.139,8.314,9.463,
*    &		             10.614,13.034,15.434,17.859/
*	real*4    engch2(16)/1.674,2.249,2.860,3.479,
*    &		             4.100,4.750,5.338,5.920,
*    &		             7.139,8.314,9.463,10.614,
ce---&		             13.034,15.434,17.859,20.284/
	real*4	  sdata(nlim,2),esdata(nlim,2)
	real*4	  sigma(2),esigma(2)
        real*4    ratio(10,nlim),rerr(10,nlim)
	real*8	  mjdhjm,mjdowr
	integer   ntime0(7),ntime1(7)
	integer	  nch1(10,2),nch2(10,2)
	integer   ipage
 	character cmnt*20,text(30)*40
ce
	real*4	  engch1(16),engch2(16)
        common /engrng/ engch1,engch2
*------- Clear data ------------
 1101	do 1100 i11=1,nlim
 1301	  do 1300 i13=1,10
	    ratio(i13,i11)=0.0
	    rerr(i13,i11)=0.0
 1300	  continue
 1100	continue
*------- Program Start ----------------
	num=0
 9000	continue
	num=num+1
	if(num.gt.10) goto 999 
*------- Initial ---------------
	xmax=1800.0
	ymax=0.0
	xmin=0.0
	ymin=0.0
*------- Input channel ---------
	write(*,*) '*------- Hardness (Softness) Ratio -------*'
	write(*,*)
  940	continue
*------- Input A ---------------
  941	continue
	write(*,*) 'Input channel (A1-A2/B1-B2)'
	write(*,*)
	write(*,*) 'A1 ?'
	read(*,'(i2)') nch1(num,1)
	write(*,*) 'A2 ?'
	read(*,'(i2)') nch2(num,1)
	write(*,*)
*------- Judgement -------------
	if(nch1(num,1).gt.16.or.nch2(num,1).gt.16) goto 941
	if(nch1(num,1).gt.nch2(num,1)) goto 941
*------- Input B ---------------
  942	continue
	write(*,*) 'B1 ?'
	read(*,'(i2)') nch1(num,2)
	write(*,*) 'B2 ?'
	read(*,'(i2)') nch2(num,2)
	write(*,*)
*------- Judgement -------------
	if(nch1(num,2).gt.16.or.nch2(num,2).gt.16) goto 942
	if(nch1(num,2).gt.nch2(num,2)) goto 942
*------- plot channel ----------
	write(text(num),101)
     &	engch1(nch1(num,1)),engch2(nch2(num,1)),
     &	engch1(nch1(num,2)),engch2(nch2(num,2))
  101	format(f5.2,'-',f5.2,'/',f5.2,'-',f5.2,' keV')
*------- Calculate Ratio ------
 1201	do 1200 i12=1,ndata
 1211	  do 1210 i22=1,2
	    lch1=nch1(num,i22)
	    lch2=nch2(num,i22)
	    if(lch1.eq.1.and.lch2.eq.8) then
	      lch1=17
	      lch2=17
	    elseif(lch1.eq.9.and.lch2.eq.16) then
	      lch1=18
	      lch2=18
	    elseif(lch1.eq.1.and.lch2.eq.16) then
	      lch1=19
	      lch2=19
	    endif
	    sdata(i12,i22)=0.0
	    esdata(i12,i22)=0.0
	    if(lch1.eq.lch2) then
	      sdata(i12,i22)=data(i12,lch1)
	      esdata(i12,i22)=edata(i12,lch1)
	    else
 1221	      do 1220 j21=lch1,lch2
	        if(edata(i12,j21).le.0.0) then
	          sdata(i12,i22)=-9.
	          esdata(i12,i22)=-9.
		  goto 1210
		else
	          sdata(i12,i22)=sdata(i12,i22)+data(i12,j21)
	          esdata(i12,i22)=esdata(i12,i22)**2
     &	                               +edata(i12,j21)**2
		endif
 1220	      continue
	      if(esdata(i12,i22).ge.0.0)
     &	         esdata(i12,i22)=sqrt(esdata(i12,i22))
	    endif
 1210	  continue
	  if(sdata(i12,1).le.0.0.or.sdata(i12,2).le.0.0) then
	    ratio(num,i12)=-99.
	    rerr(num,i12)=-99.
	  else
	    ratio(num,i12)=sdata(i12,1)/sdata(i12,2)
	    rerr(num,i12)=sqrt((sdata(i12,1)*esdata(i12,2))
     &	     **2+(sdata(i12,2)*esdata(i12,1))**2)/(sdata(i12,2)**2)
	  endif
 1200	continue
*------- Clear data -----------
	 sigma(1)=0.0
	 sigma(2)=0.0
	esigma(1)=0.0
	esigma(2)=0.0
	 xmean=0.0
	exmean=0.0
*------- Calcurate Mean -----------
 1401	do 1400 i14=1,ndata
	  if(esdata(i14,1).le.0.0.or.esdata(i14,2).le.0.0)
     &	     goto 1400
	  sigma(1)=sigma(1)+sdata(i14,1)/(esdata(i14,1)**2)
	  sigma(2)=sigma(2)+sdata(i14,2)/(esdata(i14,2)**2)
	  esigma(1)=esigma(1)+1.0/(esdata(i14,1)**2)
	  esigma(2)=esigma(2)+1.0/(esdata(i14,2)**2)
 1400	continue
*------- Mean value ---------------
	if(sigma(1).gt.0.0.and.sigma(2).gt.0.0) then
	  xmean=(sigma(1)*esigma(2))/(sigma(2)*esigma(1))
	 exmean=(sqrt((sigma(2)**2)*esigma(1)+(sigma(1)**2)
     &              *esigma(2))*esigma(2))/((sigma(2)**2)*esigma(1))
	  write(*,129) xmean,exmean
  129	  format(' Mean value = ',F8.4,' +-',F8.4)
	  write(*,*)
	else
	  write(*,*) 'Mean value cannot be calcurated...'
	  write(*,*)
	endif
*------- Scale --------------------
	write(*,*) 'Plot all you inputted (y-1/n-0) ?'
	read(*,'(i2)') id
	write(*,*)
*------- Y - jiku -----------------
	ymax1=ratio(num,1)
	if(id.eq.1) then
 5001	  do 5000 i51=1,num
 5101	    do 5100 i52=1,ndata
	      if(ymax1.lt.ratio(i51,i52)) ymax1=ratio(i51,i52)
 5100	    continue
 5000	  continue
	else
 5201	  do 5200 i53=1,ndata
	    if(ymax1.lt.ratio(num,i53)) ymax1=ratio(num,i53)
 5200	  continue
	endif
	if(ymax1.le.0.0) ymax1=1.0
*------ Waku wo kaku ------------------------
        if(ymax1.lt.1.0) then
          ymax1=ymax1*100000.
          iflag=1
        endif
        call sutble(ymax1,ymax)
        if(iflag.eq.1) ymax=ymax/100000.
	xmin=0.0
	ymin=0.0
	xmax0=xmax
	ymax0=ymax
*------- Plot -------------------------------
 8800	continue
*------- Calcurate MJD range ------
ce	mjdhjm=dble(xmin)+46795.5d0
ce	mjdowr=dble(xmax)+46795.5d0
	mjdhjm=dble(xmin)+46796.0d0
	mjdowr=dble(xmax)+46796.0d0
	call mjdate(mjdhjm,ntime0)
	call mjdate(mjdowr,ntime1)
*------- Make Comment -------------
	write(text(21),600) (ntime0(j),j=1,3),(ntime1(k),k=1,3)
  600	format('  from ',i4,2('/ ',i2),'  to ',i4,2('/ ',i2))
	call rplot
     &  (num,day,data,edata,ndata,ratio,rerr,
     &	 text,cmnt,xmin,xmax,ymin,ymax,id,ipage)
*------- Change Scale -------------
	write(*,*)
	write(*,*) 'Change the scales (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 8900
	call scalet(xmin,xmax,ymin,ymax)
	if(xmax.le.0.0) xmax=xmax0
	if(ymax.le.0.0) ymax=ymax0	
	goto 8800
*------- Next -------------------------------
 8900	continue
	write(*,*) 'Plot more (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.eq.1) goto 9000
*------ Return -----------------
  999	continue
	return
	end
********************************************************
	subroutine rplot
     &  (num,day,data,edata,ndata,ratio,rerr,
     &	 text,cmnt,xmin,xmax,ymin,ymax,id,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	real*4	  ratio(10,nlim),rerr(10,nlim)
 	character cmnt*20,text(30)*40
	integer   ipage
*------ Plot Start -------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
*------ Waku -------------------
	call pgsch(0.8)
	call pgvport(0.1,0.9,0.52,0.92)
	call pgscr(19,0.8,0.6,0.1)
	call pgsci(19)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin,ymax)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(17,0.2,1.0,0.4)
	call pgsci(17)
	call pgscf(2)
	call pgsch(0.95)
	call pglabel(' ','Ratio','Hardness Ratio of '//cmnt(1:20))
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.8)
	ispec=1
	if(id.eq.1) ispec=num
 2001	do 2000 i2=1,ispec
	  ispp=i2
	  if(id.ne.1) ispp=num
 2201	  do 2200 i21=1,ndata
	    if(ratio(ispp,i21).le.0.0) goto 2200
	    x0=day(i21)
	    y0=ratio(ispp,i21)
	    y1=ratio(ispp,i21)-rerr(ispp,i21)
	    y2=ratio(ispp,i21)+rerr(ispp,i21)
	    call pgsci(ispp+1)
	    call pgpoint(1,x0,y0,ispp+1)
	    call pgerry(1,x0,y1,y2,0.0)
 2200	  continue
 2000	continue
*------ Light Curve (1-20 keV) -
	call pgvport(0.1,0.9,0.22,0.47)
	call pgsch(0.8)
*------ Y jiku no scale --------
	ymin2=0.0
	ymax2=data(1,19)
 3001	do 3000 i31=2,ndata
	  if(ymax2.lt.data(i31,19)) ymax2=data(i31,19)
 3000	continue
*------ Waku wo kaku -----------
	iflag=1
	if(ymax2.le.1.0) then
	  ymax2=ymax2*100000.
	  iflag=4
	endif
	call sutble(ymax2,ymax3)
	if(iflag.eq.4) ymax3=ymax3/100000.
	call pgsci(19)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin2,ymax3)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(18,1.0,0.6,0.2)
	call pgsci(17)
	call pgscf(2)
	call pgsch(0.95)
	call pglabel('Day of 1987','Flux',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.7)
 4001	do 4000 i41=1,ndata
	  if(data(i41,19).le.0.0) goto 4000
	  x6=day(i41)
	  y6=data(i41,19)
	  y7=data(i41,19)-edata(i41,19)
	  y8=data(i41,19)+edata(i41,19)
	  call pgsci(7)
	  call pgpoint(1,x6,y6,5)
	  call pgerry(1,x6,y7,y8,0.0)
 4000	continue
*------ title ------------------
	call pgscr(21,1.0,0.6,0.7)
	call pgsci(21)
	call pgscf(3)
	call pgsch(0.8)
	x9=xmin+(xmax-xmin)*0.85
	y9=ymin2+(ymax3-ymin2)*0.8
	call pgtext(x9,y9,'1-20 keV')
*------ Plot coment ------------
	call pgvport(0.1,0.9,0.0,0.12)
	call pgwindow(0.0,1.0,0.0,1.0)
	call pgsch(0.75)
	call pgscf(3)
	x4=0.02
	y4=0.9
	call pgscr(20,0.5,1.0,1.0)
	call pgsci(20)
	call pgtext(x4,y4,text(21))
	call pgscf(1)
	ispr=1
	if(id.eq.1) ispr=num
 7001	do 7000 i7=1,ispr
	  ispt=i7
	  if(id.ne.1) ispt=num
	  call pgsci(ispt+1)
	  y4=y4-0.22
	  if(y4.lt.0.05) then
	    x4=x4+0.32
	    if(x4.gt.0.9) goto 7000
	    y4=0.68
	  endif
	  call pgtext(x4,y4,char(ispt+1)//text(ispt))
 7000	continue
*------- Return -----------------
  999	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
********************************************************
	subroutine lcurve(day,data,edata,ndata,cmnt,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
ce---	real*4    engch1(16)/1.000,1.674,2.249,2.860,
*    &		             3.479,4.100,4.750,5.338,
*    &		             5.920,7.139,8.314,9.463,
*    &		             10.614,13.034,15.434,17.859/
*	real*4    engch2(16)/1.674,2.249,2.860,3.479,
*    &		             4.100,4.750,5.338,5.920,
*    &		             7.139,8.314,9.463,10.614,
ce---&		             13.034,15.434,17.859,20.284/
	real*4	  sdata(10,nlim),esdata(10,nlim)
	real*8	  mjdhjm,mjdowr
	integer   ntime0(7),ntime1(7)
	integer	  nch1(10),nch2(10)
	integer   ipage
 	character cmnt*20,text(30)*40
ce
	real*4	  engch1(16),engch2(16)
        common /engrng/ engch1,engch2
*------- Clear data ------------
 1101	do 1100 i11=1,nlim
 1301	  do 1300 i13=1,10
	    sdata(i13,i11)=0.0
	    esdata(i13,i11)=0.0
 1300	  continue
 1100	continue
*------- Program Start ----------------
	num=0
 9000	continue
	num=num+1
	if(num.gt.10) goto 999 
*------- Initial ---------------
	xmax=1800.0
	ymax=0.0
	xmin=0.0
	ymin=0.0
*------- Input channel ---------
	write(*,*) '*--------- ASM Light Curve ---------*'
	write(*,*)
  940	continue
*------- Input A ---------------
  941	continue
	write(*,*) 'Input channel (A1-A2)'
	write(*,*)
	write(*,*) 'A1 ?'
*	read(*,'(i2)') nch1(num)
	read(*,'(i2)') nch1(num)
	write(*,*) 'A2 ?'
*	read(*,'(i2)') nch2(num)
	read(*,'(i2)') nch2(num)
	write(*,*)
*------- Judgement -------------
	if(nch1(num).gt.16.or.nch2(num).gt.16) goto 941
	if(nch1(num).gt.nch2(num)) goto 941
*------- plot channel ----------
	write(text(num),101) engch1(nch1(num)),engch2(nch2(num))
  101	format(f5.2,'-',f5.2,' keV')
*------- Calculate data --------
	xdata=0.0
	exdata=0.0
 1201	do 1200 i12=1,ndata
	  lch1=nch1(num)
	  lch2=nch2(num)
	  sdata(num,i12)=0.0
	  esdata(num,i12)=0.0
	  if(lch1.eq.1.and.lch2.eq.8) then
	    lch1=17
	    lch2=17
	  elseif(lch1.eq.9.and.lch2.eq.16) then
	    lch1=18
	    lch2=18
	  elseif(lch1.eq.1.and.lch2.eq.16) then
	    lch1=19
	    lch2=19
	  endif
	  if(lch1.eq.lch2) then
	    sdata(num,i12)=data(i12,lch1)
	    esdata(num,i12)=edata(i12,lch1)
	  else
 1221	    do 1220 j21=lch1,lch2
	      if(edata(i12,j21).le.0.0) then
	        sdata(num,i12)=-9.
	        esdata(num,i12)=-9.
		goto 1200
	      else
	        sdata(num,i12)=sdata(num,i12)+data(i12,j21)
	        esdata(num,i12)=esdata(num,i12)**2
     &	                             +edata(i12,j21)**2
	      endif
 1220	    continue
	    if(esdata(num,i12).ge.0.0)
     &	       esdata(num,i12)=sqrt(esdata(num,i12))
	  endif
	  if(esdata(num,i12).gt.0.0) then
	    xdata=xdata+sdata(num,i12)/(esdata(num,i12)**2)
	    exdata=exdata+1.0/(esdata(num,i12)**2)
	  endif
 1200	continue
*------- Mean value ---------------
	xdata=xdata/exdata
	exdata=1.0/sqrt(exdata)
	write(*,129) xdata,exdata
  129	format(' Mean value = ',F8.4,' +-',F8.4)
	write(*,*)
*------- Scale --------------------
	write(*,*) 'Plot all you inputted (y-1/n-0) ?'
	read(*,'(i2)') id
	write(*,*)
*------- Y - jiku -----------------
	ymax1=sdata(num,1)
	if(id.eq.1) then
 5001	  do 5000 i51=1,num
 5101	    do 5100 i52=1,ndata
	      if(ymax1.lt.sdata(i51,i52)) ymax1=sdata(i51,i52)
 5100	    continue
 5000	  continue
	else
 5201	  do 5200 i53=1,ndata
	    if(ymax1.lt.sdata(num,i53)) ymax1=sdata(num,i53)
 5200	  continue
	endif
	if(ymax1.le.0.0) ymax1=1.0
*------ Waku wo kaku ------------------------
        if(ymax1.lt.1.0) then
          ymax1=ymax1*100000.
          iflag=1
        endif
        call sutble(ymax1,ymax)
        if(iflag.eq.1) ymax=ymax/100000.
	xmin=0.0
	ymin=0.0
	xmax0=xmax
	ymax0=ymax
*------- Plot -------------------------------
 8800	continue
*------- Calcurate MJD range ------
ce	mjdhjm=dble(xmin)+46795.5d0
ce	mjdowr=dble(xmax)+46795.5d0
	mjdhjm=dble(xmin)+46796.0d0
	mjdowr=dble(xmax)+46796.0d0
	call mjdate(mjdhjm,ntime0)
	call mjdate(mjdowr,ntime1)
*------- Make Comment -------------
	write(text(21),600) (ntime0(j),j=1,3),(ntime1(k),k=1,3)
  600	format('  from ',i4,2('/ ',i2),'  to ',i4,2('/ ',i2))
	write(*,*) 'Linear (1) or Log (2) ?'
	read(*,'(i2)') ill
	write(*,*)
 	if(ill.eq.1) then
	  call lplot1
     &    (num,day,data,edata,ndata,sdata,esdata,
     &	   text,cmnt,xmin,xmax,ymin,ymax,id,ipage)
	else
	  call lplot2
     &    (num,day,data,edata,ndata,sdata,esdata,
     &	   text,cmnt,xmin,xmax,ymin,ymax,id,ipage)
	endif
*------- Change Scale -------------
	write(*,*)
	write(*,*) 'Change the scales (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 8900
	call scalet(xmin,xmax,ymin,ymax)
	if(xmax.le.0.0) xmax=xmax0
	if(ymax.le.0.0) ymax=ymax0	
	goto 8800
*------- Next -------------------------------
 8900	continue
	write(*,*) 'Plot more (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.eq.1) goto 9000
*------ Return -----------------
  999	continue
	return
	end
********************************************************
	subroutine lplot1
     &  (num,day,data,edata,ndata,sdata,esdata,
     &	 text,cmnt,xmin,xmax,ymin,ymax,id,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	real*4	  sdata(10,nlim),esdata(10,nlim)
 	character cmnt*20,text(30)*40
	integer   ipage
*------ Plot Start -------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
*------ Waku -------------------
	call pgsch(0.9)
	call pgvport(0.05,0.97,0.1,0.68)
	call pgscr(19,0.8,0.8,0.8)
	call pgsci(19)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin,ymax)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgsch(0.93)
	call pgscr(17,0.2,1.0,0.4)
	call pgsci(17)
	call pgscf(2)
ce	call pglabel('Day of 1987','Flux (counts/sec/cm**2)',' ')
	call pglabel('Day of 1987','Flux (counts/sec/cm\\u2\\d)',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.83)
	ispec=1
	if(id.eq.1) ispec=num
 2001	do 2000 i2=1,ispec
	  ispp=i2
	  if(id.ne.1) ispp=num
 2201	  do 2200 i21=1,ndata
	    if(esdata(ispp,i21).le.0.0) goto 2200
	    x0=day(i21)
	    y0=sdata(ispp,i21)
	    y1=sdata(ispp,i21)-esdata(ispp,i21)
	    y2=sdata(ispp,i21)+esdata(ispp,i21)
	    call pgsci(ispp+1)
	    call pgpoint(1,x0,y0,ispp+1)
	    call pgerry(1,x0,y1,y2,0.0)
 2200	  continue
 2000	continue
*------ Light Curve (1-20 keV) -
	call pgvport(0.55,0.97,0.78,0.97)
	call pgsch(0.75)
*------ Y jiku no scale --------
	ymin2=0.0
	ymax2=data(1,19)
 3001	do 3000 i31=2,ndata
	  if(ymax2.lt.data(i31,19)) ymax2=data(i31,19)
 3000	continue
*------ Waku wo kaku -----------
	iflag=1
	if(ymax2.le.1.0) then
	  ymax2=ymax2*100000.
	  iflag=4
	endif
	call sutble(ymax2,ymax3)
	if(iflag.eq.4) ymax3=ymax3/100000.
	call pgsci(19)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin2,ymax3)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(18,1.0,0.6,0.2)
	call pgsci(18)
	call pgscf(2)
	call pgsch(0.78)
	call pglabel('Day of 1987','Flux',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.6)
 4001	do 4000 i41=1,ndata
	  if(data(i41,19).le.0.0) goto 4000
	  x6=day(i41)
	  y6=data(i41,19)
	  y7=data(i41,19)-edata(i41,19)
	  y8=data(i41,19)+edata(i41,19)
	  call pgsci(7)
	  call pgpoint(1,x6,y6,7)
	  call pgerry(1,x6,y7,y8,0.0)
 4000	continue
*------ title ------------------
	call pgscr(21,1.0,0.6,0.7)
	call pgsci(21)
	call pgscf(3)
	call pgsch(0.7)
	x9=xmin+(xmax-xmin)*0.8
	y9=ymin2+(ymax3-ymin2)*0.8
	call pgtext(x9,y9,'1-20 keV')
*------ Plot coment ------------
	call pgvport(0.05,0.5,0.72,0.97)
	call pgwindow(0.0,1.0,0.0,1.0)
	call pgsci(13)
	call pgsch(0.97)
	call pgscf(2)
	call pgtext(0.05,0.95,'ASM Light Curve of '//cmnt(1:20))
	call pgsch(0.8)
	call pgscf(3)
	x4=0.02
	y4=0.75
	call pgscr(20,0.5,1.0,1.0)
	call pgsci(20)
	call pgtext(x4,y4,text(21))
	call pgscf(1)
	ispr=1
	if(id.eq.1) ispr=num
 7001	do 7000 i7=1,ispr
	  ispt=i7
	  if(id.ne.1) ispt=num
	  call pgsci(ispt+1)
	  y4=y4-0.14
	  if(y4.lt.0.05) then
	    x4=0.52
	    y4=0.61
	  endif
	  call pgtext(x4,y4,char(ispt+1)//text(ispt)(1:17))
 7000	continue
*------- Return -----------------
  999	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
********************************************************
	subroutine lplot2
     &  (num,day,data,edata,ndata,sdata,esdata,
     &	 text,cmnt,xmin,xmax,ymin,ymax,id,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	real*4	  sdata(10,nlim),esdata(10,nlim)
	real*8	  mjdhjm,mjdowr,mjday
 	character cmnt*20,text(30)*40
	integer   ipage
*------- Kirino ii suuji --------------------
	if(ymax.le.0.0) ymax=1.0
	iloy=int(alog10(ymax))+1
	ylog1=float(iloy-1)-alog10(ymax)
	if(ylog1.gt.0.0) iloy=iloy-1
	ymax1=iloy
	if(ymin.le.0.0) then
	  ymin1=iloy-3
	  ymin0=10**ymin1
	else
	  ymin0=ymin
	  ymin1=int(alog10(ymin))
	endif
*------ Plot Start -------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
*------ Waku -------------------
	call pgsch(0.9)
	call pgvport(0.05,0.97,0.1,0.68)
	call pgscr(19,0.8,0.8,0.8)
	call pgsci(19)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin1,ymax1)
	call pgbox('BCTSN',0.0,0,'LBCTSN',0.0,0)
*------ plot comment -----------
	call pgsch(0.93)
	call pgscr(17,0.2,1.0,0.4)
	call pgsci(17)
	call pgscf(2)
ce	call pglabel('Day of 1987','Flux (counts/sec/cm**2)',' ')
	call pglabel('Day of 1987','Flux (counts/sec/cm\\u2\\d)',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.83)
	ispec=1
	if(id.eq.1) ispec=num
 2001	do 2000 i2=1,ispec
	  ispp=i2
	  if(id.ne.1) ispp=num
 2201	  do 2200 i21=1,ndata
	    x0=day(i21)
	    if(sdata(ispp,i21).le.0.0) then
	      y0=alog10(ymin0*0.1)
	    else
	      y0=alog10(amax1(sdata(ispp,i21),ymin0))
	    endif
	    y1=alog10(amax1(sdata(ispp,i21)
     &		     -esdata(ispp,i21),ymin0))
	    y2=alog10(amax1(sdata(ispp,i21)
     &		     +esdata(ispp,i21),ymin0))
	    call pgsci(ispp+1)
	    call pgpoint(1,x0,y0,ispp+1)
	    call pgerry(1,x0,y1,y2,0.0)
 2200	  continue
 2000	continue
*------ Light Curve (1-20 keV) -
	call pgvport(0.55,0.97,0.78,0.97)
	call pgsch(0.75)
*------ Y jiku no scale --------
	ymin2=0.0
	ymax2=data(1,19)
 3001	do 3000 i31=2,ndata
	  if(ymax2.lt.data(i31,19)) ymax2=data(i31,19)
 3000	continue
*------- Kirino ii suuji --------------------
	iloy=int(alog10(ymax2))+1
	ylog2=float(iloy-1)-alog10(ymax2)
	if(ylog2.gt.0.0) iloy=iloy-1
	ymin2=iloy-3
	ymin0=10**ymin2
	ymax3=iloy
	call pgsci(19)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin2,ymax3)
	call pgbox('BCTSN',0.0,0,'LBCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(18,1.0,0.6,0.2)
	call pgsci(18)
	call pgscf(2)
	call pgsch(0.78)
	call pglabel('Day of 1987','Flux',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.6)
 4001	do 4000 i41=1,ndata
	  x6=day(i41)
	  if(data(i41,19).le.0.0) then
	    y6=alog10(ymin0*0.1)
	  else
	    y6=alog10(amax1(data(i41,19),ymin0))
	  endif
	  y7=alog10(amax1(data(i41,19)-edata(i41,19),ymin0))
	  y8=alog10(amax1(data(i41,19)+edata(i41,19),ymin0))
	  call pgsci(7)
	  call pgpoint(1,x6,y6,7)
	  call pgerry(1,x6,y7,y8,0.0)
 4000	continue
*------ title ------------------
	call pgscr(21,1.0,0.6,0.7)
	call pgsci(21)
	call pgscf(3)
	call pgsch(0.7)
	x9=xmin+(xmax-xmin)*0.8
	y9=ymin2+(ymax3-ymin2)*0.8
	call pgtext(x9,y9,'1-20 keV')
*------ Plot coment ------------
	call pgvport(0.05,0.5,0.72,0.97)
	call pgwindow(0.0,1.0,0.0,1.0)
	call pgsci(13)
	call pgsch(0.97)
	call pgscf(2)
	call pgtext(0.05,0.95,'ASM Light Curve of '//cmnt(1:20))
	call pgsch(0.8)
	call pgscf(3)
	x4=0.02
	y4=0.75
	call pgscr(20,0.5,1.0,1.0)
	call pgsci(20)
	call pgtext(x4,y4,text(21))
	call pgscf(1)
	ispr=1
	if(id.eq.1) ispr=num
 7001	do 7000 i7=1,ispr
	  ispt=i7
	  if(id.ne.1) ispt=num
	  call pgsci(ispt+1)
	  y4=y4-0.14
	  if(y4.lt.0.05) then
	    x4=0.52
	    y4=0.61
	  endif
	  call pgtext(x4,y4,char(ispt+1)//text(ispt)(1:17))
 7000	continue
*------- Return -----------------
  999	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
****************************************************************
        subroutine color(nch,day,data,edata,ndata,cmnt,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
        real*4    r1(10,nlim),re1(10,nlim)
	real*4	  r2(10,nlim),re2(10,nlim)
ce---	real*4    engch1(16)/1.000,1.674,2.249,2.860,
*    &		             3.479,4.100,4.750,5.338,
*    &		             5.920,7.139,8.314,9.463,
*    &		             10.614,13.034,15.434,17.859/
*	real*4    engch2(16)/1.674,2.249,2.860,3.479,
*    &		             4.100,4.750,5.338,5.920,
*    &		             7.139,8.314,9.463,10.614,
ce---&		             13.034,15.434,17.859,20.284/
	real*8	  mjdhjm,mjdowr,mjday
	real*4	  sigma(4),esigma(4)
	real*4	  sdata(nlim,4),esdata(nlim,4)
	real*4	  sttmjd(10),stpmjd(10)
	real*4	  lday(10,nlim),ldata(10,nlim),ledata(10,nlim)
	integer	  lndata(10)
	integer   nch(8)
	integer	  id,ipage
	dimension ntime0(7),ntime1(7)
 	character cmnt*20,text(30)*50
	character text1*47,text2*47
ce
	real*4	  engch1(16),engch2(16)
        common /engrng/ engch1,engch2
*------- plot channel ----------
	write(text1,101)
     &  engch1(nch(1)),engch2(nch(2)),engch1(nch(3)),engch2(nch(4))
	write(text2,101)
     &  engch1(nch(5)),engch2(nch(6)),engch1(nch(7)),engch2(nch(8))
  101	format(f5.2,'-',f5.2,'/',f5.2,'-',f5.2,' keV')
*------- Clear data ------------
 1101	do 1100 i11=1,nlim
 1301	  do 1300 i13=1,10
	    r1(i13,i11)=0.0
	    r2(i13,i11)=0.0
	    re1(i13,i11)=0.0
	    re2(i13,i11)=0.0
 1300	  continue
 1100	continue
*------- Program Start ----------------
	num=0
ce------ Data Selection Mode ----------
 200    continue
        write(*,*) 'Please input data selection mode'
        write(*,*) '1...Day  2...Flux  3...Phase'
        read(*,'(i2)') imode
        write(*,*)
        if (imode.lt.1 .or. imode.gt.3) goto 200
        if (imode.eq.3) then
           write(*,*) 'Input Period (Day)'
           read(*,*) period
           write(*,*)
        endif
        mjdhjm=46796.0d0
        mjdowr=46796.0d0+1800.0d0
ce---
 9000	continue
	num=num+1
	if(num.gt.10) goto 999 
*------- Initial ---------------
	xmax=0.0
	ymax=0.0
	xmin=0.0
	ymin=0.0
*------- Read day to add data ---------
ce	write(*,*) 'Input day from 1987 Jan.1'
*	write(*,*) 
*	write(*,*) 'Start?'
*	read(*,*) mjdhjm
*ce	mjdhjm=mjdhjm+46795.5d0
*	mjdhjm=mjdhjm+46796.0d0
*	write(*,*) 
*	write(*,*) 'Stop?'
*	read(*,*) mjdowr
*ce	mjdowr=mjdowr+46795.5d0
*	mjdowr=mjdowr+46796.0d0
ce	write(*,*) 
ce------ Read Day or Discrimination level to add data -------
        if (imode.eq.1) then
 1001      do 1000 i1=4,7
              ntime0(i1)=0
              ntime1(i1)=0
 1000      continue
           write(*,*) 'Input Start day.'
           write(*,*)
           write(*,*) 'Year or Day from 1987 Jan 1 ?'
           read(*,*) sday
           if (sday.lt.1987 .or. sday.gt.1992) then
              mjdhjm=dble(sday)+46796.0d0
              write(*,*)
              call mjdate(mjdhjm,ntime0)
              goto 800
           endif
           ntime0(1)=sday
           write(*,*) 'Month ?'
*           read(*,'(i2)') ntime0(2)
           read(*,'(i2)') ntime0(2)
           write(*,*) 'Day ?'
*           read(*,*) ntime0(3)
           read(*,'(i2)') ntime0(3)
           write(*,*) 
           call mjulia(ntime0,mjdhjm)
 800       continue
           write(*,*) 'Input Stop day.'
           write(*,*)
           write(*,*) 'Year or Day from 1987 Jan 1 ?'
           read(*,*) sday
           if (sday.lt.1987 .or. sday.gt.1992) then
              mjdowr=dble(sday)+46796.0d0
              write(*,*)
              call mjdate(mjdowr,ntime1)
              goto 1400
           endif
           ntime1(1)=sday
           write(*,*) 'Month ?'
*           read(*,*) ntime1(2)
           read(*,'(i2)') ntime1(2)
           write(*,*) 'Day ?'
*           read(*,*) ntime1(3)
           read(*,'(i2)') ntime1(3)
           write(*,*) 
           call mjulia(ntime1,mjdowr)
        elseif (imode.eq.2) then
           write(*,*) 'Input Upper Discrimination.'
           read(*,*) discu
           write(*,*)
           write(*,*) 'Input Lower Discrimination.'
           read(*,*) discl
           write(*,*)
           call mjdate(mjdhjm,ntime0)
           call mjdate(mjdowr,ntime1)
        elseif (imode.eq.3) then
           write(*,*) 'Input Start Phase.'
           read(*,*) stpha
           write(*,*)
           write(*,*) 'Input Stop Phase.'
           read(*,*) sppha
           write(*,*)
           call mjdate(mjdhjm,ntime0)
           call mjdate(mjdowr,ntime1)
        endif
*------- Data no hajime to Owari -------
ce	call mjdate(mjdhjm,ntime0)
ce	call mjdate(mjdowr,ntime1)
 1400 continue
	sttmjd(num)=sngl(mjdhjm)
	stpmjd(num)=sngl(mjdowr)
*------- Clear data -----------
 1501	do 1500 i15=1,4
	  sigma(i15)=0.0
	  esigma(i15)=0.0
 1500	continue
	xm1=0.0
	exm1=0.0
	xm2=0.0
	exm2=0.0
	mndata=0
	lndata(num)=0
*------- Calculate Ratio ------
 1201	do 1200 i12=1,ndata
ce	  mjday=dble(day(i12))+46795.5d0
	  mjday=dble(day(i12))+46796.0d0
ce---	  if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 1200
          if(imode.eq.1) then
             if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 1200
          elseif(imode.eq.2) then
             if(data(i12,17).lt.discl .or. data(i12,17).gt.discu)
     &            goto 1200
          elseif(imode.eq.3) then
             phase = day(i12)/period - int(day(i12)/period)
             if(sppha .gt. stpha) then
                if(phase.lt.stpha .or. phase.gt.sppha)
     &               goto 1200
             elseif(sppha .lt. stpha) then
                if(phase.lt.stpha .and. phase.gt.sppha)
     &               goto 1200
             endif
          endif
ce---
	  mndata=mndata+1
	  lday(num,mndata)=sngl(mjday)
	  ldata(num,mndata)=data(i12,19)
	  ledata(num,mndata)=edata(i12,19)
 1211	  do 1210 i22=1,4
	    lch1=nch(2*i22-1)
	    lch2=nch(2*i22)
	    sdata(mndata,i22)=0.0
	    esdata(mndata,i22)=0.0
	    if(lch1.eq.1.and.lch2.eq.8) then
		lch1=17
		lch2=17
	    elseif(lch1.eq.9.and.lch2.eq.16) then
		lch1=18
		lch2=18
	    elseif(lch1.eq.1.and.lch2.eq.16) then
		lch1=19
		lch2=19
	    endif
 1221	    do 1220 j21=lch1,lch2
	      if(edata(i12,j21).le.0.0) then
	        sdata(mndata,i22)=-9.
	        esdata(mndata,i22)=-9.
		goto 1210
	      else
	        sdata(mndata,i22)=sdata(mndata,i22)+data(i12,j21)
	        esdata(mndata,i22)=esdata(mndata,i22)**2
     &	                                   +edata(i12,j21)**2
	      endif
 1220	    continue
	    if(esdata(mndata,i22).gt.0.0)
     &	       esdata(mndata,i22)=sqrt(esdata(mndata,i22))
 1210	  continue
	  if(esdata(mndata,1).gt.0.0.and.
     &      esdata(mndata,2).gt.0.0) then
	    sigma(1)=sigma(1)+sdata(mndata,1)/(esdata(mndata,1)**2)
	    esigma(1)=esigma(1)+1.0/(esdata(mndata,1)**2)
	    sigma(2)=sigma(2)+sdata(mndata,2)/(esdata(mndata,2)**2)
	    esigma(2)=esigma(2)+1.0/(esdata(mndata,2)**2)
	  endif
	  if(esdata(mndata,3).gt.0.0.and.
     &      esdata(mndata,4).gt.0.0) then
	    sigma(3)=sigma(3)+sdata(mndata,3)/(esdata(mndata,3)**2)
	    esigma(3)=esigma(3)+1.0/(esdata(mndata,3)**2)
	    sigma(4)=sigma(4)+sdata(mndata,4)/(esdata(mndata,4)**2)
	    esigma(4)=esigma(4)+1.0/(esdata(mndata,4)**2)
	  endif
	  if(sdata(mndata,1).le.0.0.or.sdata(mndata,2).le.0.0) then
	    mndata=mndata-1
	    goto 1200
	  else
	    r1(num,mndata)=sdata(mndata,1)/sdata(mndata,2)
	    re1(num,mndata)=sqrt((sdata(mndata,1)*esdata(mndata,2))
     &	     **2+(sdata(mndata,2)*esdata(mndata,1))**2)
     &	     /(sdata(mndata,2)**2)
	  endif
	  if(sdata(mndata,3).le.0.0.or.sdata(mndata,4).le.0.0) then
	    mndata=mndata-1
	    goto 1200
	  else
	    r2(num,mndata)=sdata(mndata,3)/sdata(mndata,4)
	    re2(num,mndata)=sqrt((sdata(mndata,3)*esdata(mndata,4))
     &	     **2+(sdata(mndata,4)*esdata(mndata,3))**2)
     &	     /(sdata(mndata,4)**2)
	  endif
	  lndata(num)=mndata
 1200	continue
	if(lndata(num).le.0) then
	  write(*,*) 'No data in the period.'
	  write(*,*)
	  num=num-1
	  goto 8900
	endif
*------- Mean value ---------------
	if(sigma(1).le.0.0.and.sigma(2).le.0.0) then
	   xm1=-99.
	  exm1=-99.
	else
	   xm1=(sigma(1)*esigma(2))/(sigma(2)*esigma(1))
	  exm1=(sqrt((sigma(2)**2)*esigma(1)+(sigma(1)**2)
     &         *esigma(2))*esigma(2))/((sigma(2)**2)*esigma(1))
	endif
	if(sigma(3).le.0.0.and.sigma(4).le.0.0) then
	   xm2=-99.
	  exm2=-99.
	else
	   xm2=(sigma(3)*esigma(4))/(sigma(4)*esigma(3))
	  exm2=(sqrt((sigma(4)**2)*esigma(3)+(sigma(3)**2)
     &         *esigma(4))*esigma(4))/((sigma(4)**2)*esigma(3))
	endif
	write(*,128) xm1,exm1
  128	format(' Mean value (horizontal color) =',f8.4,' +-',f8.4)
	write(*,*)
	write(*,129) xm2,exm2
  129	format(' Mean value (vertical color) =',f8.4,' +-',f8.4)
	write(*,*)
*------- Make Comment -----------------------
	write(text(num),600) (ntime0(j),j=1,3),(ntime1(k),k=1,3)
  600	format('  From ',i4,2('/ ',i2),'  To ',i4,2('/ ',i2))
	write(text(num+10),610) lndata(num)
  610	format(8x,'The number of data points = ',i4)
*------- Scale --------------------
	write(*,*) 'Plot all you inputted (y-1/n-0) ?'
*	read(*,*) id
	read(*,'(i2)') id
	write(*,*)
*------- X - jiku -----------------
	xmax1=r1(num,1)
	if(id.eq.1) then
 4001	  do 4000 i41=1,num
 4101	    do 4100 i42=1,lndata(i41)
	      if(xmax1.lt.r1(i41,i42)) xmax1=r1(i41,i42)
 4100	    continue
 4000	  continue
	else
 4201	  do 4200 i43=2,lndata(num)
	    if(xmax1.lt.r1(num,i43)) xmax1=r1(num,i43)
 4200	  continue
	endif
	if(xmax1.le.0.0) xmax1=1.0
*------- Y - jiku -----------------
	ymax1=r2(num,1)
	if(id.eq.1) then
 5001	  do 5000 i51=1,num
 5101	    do 5100 i52=1,lndata(i51)
	      if(ymax1.lt.r2(i51,i52)) ymax1=r2(i51,i52)
 5100	    continue
 5000	  continue
	else
 5201	  do 5200 i53=2,lndata(num)
	    if(ymax1.lt.r2(num,i53)) ymax1=r2(num,i53)
 5200	  continue
	endif
	if(ymax1.le.0.0) ymax1=1.0
*------ Waku wo kaku ------------------------
        if(xmax1.lt.1.0) then
          xmax1=xmax1*100000.
          iflag=1
        endif
        call sutble(xmax1,xmax)
        if(iflag.eq.1) xmax=xmax/100000.
        if(ymax1.lt.1.0) then
          ymax1=ymax1*100000.
          iflag=2
        endif
        call sutble(ymax1,ymax)
        if(iflag.eq.2) ymax=ymax/100000.
	xmin=0.0
	ymin=0.0
	xmax0=xmax
	ymax0=ymax
*------- Plot -------------------------------
 8800	continue
	call cplot
     &  (num,r1,re1,r2,re2,text1,text2,text,cmnt,xmin,xmax,
     &	ymin,ymax,sttmjd,stpmjd,lndata,lday,ldata,ledata,
     &	id,ipage)
*------- Change Scale -------------
	write(*,*)
	write(*,*) 'Change the scales (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 8900
	call scalet(xmin,xmax,ymin,ymax)
	if(xmax.le.0.0) xmax=xmax0
	if(ymax.le.0.0) ymax=ymax0	
	goto 8800
*------- Next -------------------------------
 8900	continue
	write(*,*) 'Plot more (y-1/n-0) ?'
*	read(*,*) iyn
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.eq.1) goto 9000
*------ Return -----------------
  999	continue
	return
	end
******************************************************************
	subroutine cplot
     &  (num,r1,re1,r2,re2,text1,text2,text,cmnt,xmin,xmax,
     &	ymin,ymax,sttmjd,stpmjd,lndata,lday,ldata,ledata,
     &	id,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
        real*4    r1(10,nlim),re1(10,nlim)
	real*4	  r2(10,nlim),re2(10,nlim)
	real*4	  lday(10,200),ldata(10,200),ledata(10,200)
	real*4	  sttmjd(10),stpmjd(10)
	integer	  lndata(10)
	integer	  id,ipage
 	character cmnt*20,text(30)*50
	character text1*47,text2*47
*-------Program start -----------------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
*------- Color - Color plot -----------------
	call pgsch(0.9)
	call pgvport(0.07,0.6,0.1,0.9)
	call pgscr(25,0.6,0.7,1.0)
	call pgsci(25)
	call pgscf(1)
        call pgwindow(xmin,xmax,ymin,ymax)
        call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------- Plot Coment ------------------------
	call pgsch(0.95)
	call pgscr(17,0.3,1.0,0.1)
	call pgsci(17)
	call pgscf(2)
        call pglabel(text1,text2,'Color - Color Diagram')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.83)
	ispec=1
	if(id.eq.1) ispec=num
 6001	do 6000 i6=1,ispec
	  ispp=i6
	  if(id.ne.1) ispp=num
 6101	  do 6100 i61=1,lndata(ispp)
            if(r1(ispp,i61).le.0.0.or.r2(ispp,i61).le.0.0)
     &		 goto 6100
            x0=r1(ispp,i61)
            y0=r2(ispp,i61)
            x1=r1(ispp,i61)-re1(ispp,i61)
            x2=r1(ispp,i61)+re1(ispp,i61)
            y1=r2(ispp,i61)-re2(ispp,i61)
            y2=r2(ispp,i61)+re2(ispp,i61)
	    call pgsci(ispp+6)
	    call pgpoint(1,x0,y0,ispp+1)
            call pgerrx(1,x1,x2,y0,0.0)
            call pgerry(1,x0,y1,y2,0.0)
 6100	  continue
 6000	continue
*------ title ------------------
	call pgscr(18,1.0,0.7,0.7)
	call pgsci(18)
	call pgscf(3)
	call pgsch(0.85)
	x3=xmin+(xmax-xmin)*0.7
	y3=ymin+(ymax-ymin)*0.9
	call pgtext(x3,y3,cmnt(1:20))
*------ Light Curve (1-20 keV) -
	call pgvport(0.68,0.97,0.65,0.9)
	call pgsch(0.78)
*------ X jiku no scale --------
	if(id.eq.1.and.num.gt.1) then
	  xmin2=sttmjd(1)-10.0
	  xmax2=stpmjd(1)+10.0
 6501	  do 6500 i65=2,num
	  if(sttmjd(i65)-10.0.lt.xmin2) xmin2=sttmjd(i65)-10.0
	  if(stpmjd(i65)+10.0.gt.xmax2) xmax2=stpmjd(i65)+10.0
 6500	  continue
	else
	  xmin2=sttmjd(num)-10.0
	  xmax2=stpmjd(num)+10.0
	endif
*------ Y jiku no scale --------
	ymin2=0.0
	ymax2=ldata(num,1)
	if(id.eq.1) then
 8001	  do 8000 i81=1,num
 8101	    do 8100 i82=1,lndata(i81)
	      if(ymax2.lt.ldata(i81,i82)) ymax2=ldata(i81,i82)
 8100	    continue
 8000	  continue
	else
 8201	  do 8200 i83=1,lndata(num)
	    if(ymax2.lt.ldata(num,i83)) ymax2=ldata(num,i83)
 8200	  continue
	endif
*------ Waku wo kaku -----------
	iflag=1
	if(ymax2.le.1.0) then
	  ymax2=ymax2*100000.
	  iflag=4
	endif
	call sutble(ymax2,ymax3)
	if(iflag.eq.4) ymax3=ymax3/100000.
	call pgsci(25)
	call pgscf(1)
	call pgwindow(xmin2,xmax2,ymin2,ymax3)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(19,1.0,0.6,0.2)
	call pgsci(19)
	call pgscf(2)
ce	call pglabel('MJD','Flux (counts/sec/cm**2)',
	call pglabel('MJD','Flux (counts/sec/cm\\u2\\d)',
     &	             'Light Curve (1-20 keV)')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.65)
	ispec=1
	if(id.eq.1) ispec=num
 8301	do 8300 i83=1,ispec
	  ispq=i83
	  if(id.ne.1) ispq=num
 8401	  do 8400 i84=1,lndata(ispq)
	    if(ldata(ispq,i84).le.0.0) goto 8400
	    x6=lday(ispq,i84)
	    y6=ldata(ispq,i84)
	    y7=ldata(ispq,i84)-ledata(ispq,i84)
	    y8=ldata(ispq,i84)+ledata(ispq,i84)
	    call pgsci(ispq+6)
	    call pgpoint(1,x6,y6,ispq+1)
	    call pgerry(1,x6,y7,y8,0.0)
 8400	  continue
 8300	continue
*------ Plot coment ------------
	call pgvport(0.65,1.0,0.0,0.6)
	call pgwindow(0.0,1.0,0.0,1.0)
	call pgsch(0.8)
	call pgscf(1)
	ispr=1
	if(id.eq.1) ispr=num
 7001	do 7000 i7=1,ispr
	  ispt=i7
	  if(id.ne.1) ispt=num
	  call pgsci(ispt+6)
	  x4=0.05
	  y4=1.0-i7*0.12
	  if(y4.lt.0.1) goto 999
	  call pgtext(x4,y4,char(ispt+1)//text(ispt))
	  call pgtext(x4,y4-0.05,text(ispt+10))
 7000	continue
*------- Return -----------------
  999	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
****************************************************************
        subroutine icolor(nch,day,data,edata,ndata,cmnt,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
        real*4    r1(10,nlim),re1(10,nlim)
	real*4	  r2(10,nlim),re2(10,nlim)
ce---	real*4    engch1(16)/1.000,1.674,2.249,2.860,
*    &		             3.479,4.100,4.750,5.338,
*    &		             5.920,7.139,8.314,9.463,
*    &		             10.614,13.034,15.434,17.859/
*	real*4    engch2(16)/1.674,2.249,2.860,3.479,
*    &		             4.100,4.750,5.338,5.920,
*    &		             7.139,8.314,9.463,10.614,
ce---&		             13.034,15.434,17.859,20.284/
	real*8	  mjdhjm,mjdowr,mjday
	real*4	  sigma(3),esigma(3)
	real*4	  sdata(nlim,4),esdata(nlim,4)
	real*4	  sttmjd(10),stpmjd(10)
	real*4	  lday(10,nlim),ldata(10,nlim),ledata(10,nlim)
	integer	  lndata(10)
	integer   nch(8)
	dimension ntime0(7),ntime1(7)
 	character cmnt*20,text(30)*50
	character text1*47,text2*47
	integer   id,ipage
ce
	real*4	  engch1(16),engch2(16)
        common /engrng/ engch1,engch2
*------- plot channel ----------
	write(text1,102)
     &  engch1(nch(1)),engch2(nch(2))
	write(text2,101)
     &  engch1(nch(3)),engch2(nch(4)),engch1(nch(5)),engch2(nch(6))
  101	format(f5.2,'-',f5.2,'/',f5.2,'-',f5.2,' keV')
ce102	format(f5.2,'-',f5.2,' keV counts/sec/cm**2')
  102	format(f5.2,'-',f5.2,' keV counts/sec/cm\\u2\\d')
*------- Clear data ------------
 1101	do 1100 i11=1,nlim
 1301	  do 1300 i13=1,10
	    r1(i13,i11)=0.0
	    r2(i13,i11)=0.0
	    re1(i13,i11)=0.0
	    re2(i13,i11)=0.0
 1300	  continue
 1100	continue
*------- Program Start ----------------
	num=0
ce------ Data Selection Mode ----------
 200    continue
        write(*,*) 'Please input data selection mode'
        write(*,*) '1...Day  3...Phase'
        read(*,'(i2)') imode
        write(*,*)
        if (imode.ne.1 .and. imode.ne.3) goto 200
        if (imode.eq.3) then
           write(*,*) 'Input Period (Day)'
           read(*,*) period
           write(*,*)
        endif
        mjdhjm=46796.0d0
        mjdowr=46796.0d0+1800.0d0
ce---
 9000	continue
	num=num+1
	if(num.gt.10) goto 999 
*------- Initial ---------------
	xmax=0.0
	ymax=0.0
	xmin=0.0
	ymin=0.0
*------- Read day to add data ---------
ce	write(*,*) 'Input day from 1987 Jan.1'
*	write(*,*) 
*	write(*,*) 'Start?'
*	read(*,*) mjdhjm
*ce	mjdhjm=mjdhjm+46795.5d0
*	mjdhjm=mjdhjm+46796.0d0
*	write(*,*) 
*	write(*,*) 'Stop?'
*	read(*,*) mjdowr
*ce	mjdowr=mjdowr+46795.5d0
*	mjdowr=mjdowr+46796.0d0
ce	write(*,*) 
ce------ Read Day or Discrimination level to add data -------
        if (imode.eq.1) then
 1001      do 1000 i1=4,7
              ntime0(i1)=0
              ntime1(i1)=0
 1000      continue
           write(*,*) 'Input Start day.'
           write(*,*)
           write(*,*) 'Year or Day from 1987 Jan 1 ?'
           read(*,*) sday
           if (sday.lt.1987 .or. sday.gt.1992) then
              mjdhjm=dble(sday)+46796.0d0
              write(*,*)
              call mjdate(mjdhjm,ntime0)
              goto 800
           endif
           ntime0(1)=sday
           write(*,*) 'Month ?'
*           read(*,*) ntime0(2)
           read(*,'(i2)') ntime0(2)
           write(*,*) 'Day ?'
*           read(*,*) ntime0(3)
           read(*,'(i2)') ntime0(3)
           write(*,*) 
           call mjulia(ntime0,mjdhjm)
 800       continue
           write(*,*) 'Input Stop day.'
           write(*,*)
           write(*,*) 'Year or Day from 1987 Jan 1 ?'
           read(*,*) sday
           if (sday.lt.1987 .or. sday.gt.1992) then
              mjdowr=dble(sday)+46796.0d0
              write(*,*)
              call mjdate(mjdowr,ntime1)
              goto 1400
           endif
           ntime1(1)=sday
           write(*,*) 'Month ?'
*           read(*,*) ntime1(2)
           read(*,'(i2)') ntime1(2)
           write(*,*) 'Day ?'
*           read(*,*) ntime1(3)
           read(*,'(i2)') ntime1(3)
           write(*,*) 
           call mjulia(ntime1,mjdowr)
        elseif (imode.eq.2) then
           write(*,*) 'Input Upper Discrimination.'
           read(*,*) discu
           write(*,*)
           write(*,*) 'Input Lower Discrimination.'
           read(*,*) discl
           write(*,*)
           call mjdate(mjdhjm,ntime0)
           call mjdate(mjdowr,ntime1)
        elseif (imode.eq.3) then
           write(*,*) 'Input Start Phase.'
           read(*,*) stpha
           write(*,*)
           write(*,*) 'Input Stop Phase.'
           read(*,*) sppha
           write(*,*)
           call mjdate(mjdhjm,ntime0)
           call mjdate(mjdowr,ntime1)
        endif
*------- Data no hajime to Owari -------
ce	call mjdate(mjdhjm,ntime0)
ce	call mjdate(mjdowr,ntime1)
 1400 continue
	sttmjd(num)=sngl(mjdhjm)
	stpmjd(num)=sngl(mjdowr)
*------- Clear data -----------
 1501	do 1500 i15=1,3
	  sigma(i15)=0.0
	  esigma(i15)=0.0
 1500	continue
	 xm1=0.0
	exm1=0.0
	 xm2=0.0
	exm2=0.0
	mndata=0
	lndata(num)=0
*------- Calculate Ratio ------
 1201	do 1200 i12=1,ndata
ce	  mjday=dble(day(i12))+46795.5d0
	  mjday=dble(day(i12))+46796.0d0
ce---	  if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 1200
          if(imode.eq.1) then
             if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 1200
          elseif(imode.eq.3) then
             phase = day(i12)/period - int(day(i12)/period)
             if(sppha .gt. stpha) then
                if(phase.lt.stpha .or. phase.gt.sppha)
     &               goto 1200
             elseif(sppha .lt. stpha) then
                if(phase.lt.stpha .and. phase.gt.sppha)
     &               goto 1200
             endif
          endif
ce---
	  mndata=mndata+1
	  lday(num,mndata)=sngl(mjday)
	  ldata(num,mndata)=data(i12,19)
	  ledata(num,mndata)=edata(i12,19)
 1211	  do 1210 i22=1,3
	    lch1=nch(2*i22-1)
	    lch2=nch(2*i22)
	    sdata(mndata,i22)=0.0
	    esdata(mndata,i22)=0.0
	    if(lch1.eq.1.and.lch2.eq.8) then
		lch1=17
		lch2=17
	    elseif(lch1.eq.9.and.lch2.eq.16) then
		lch1=18
		lch2=18
	    elseif(lch1.eq.1.and.lch2.eq.16) then
		lch1=19
		lch2=19
	    endif
 1221	    do 1220 j21=lch1,lch2
	      if(edata(i12,j21).le.0.0) then
	        sdata(mndata,i22)=-9.
	        esdata(mndata,i22)=-9.
		goto 1210
	      else
	        sdata(mndata,i22)=sdata(mndata,i22)+data(i12,j21)
	        esdata(mndata,i22)=esdata(mndata,i22)**2
     &	                                   +edata(i12,j21)**2
	      endif
 1220	    continue
	    if(esdata(mndata,i22).gt.0.0)
     &	       esdata(mndata,i22)=sqrt(esdata(mndata,i22))
 1210	  continue
	  if(esdata(mndata,1).gt.0.0) then
	    sigma(1)=sigma(1)+sdata(mndata,1)/(esdata(mndata,1)**2)
	    esigma(1)=esigma(1)+1.0/(esdata(mndata,1)**2)
	  endif
	  if(esdata(mndata,2).gt.0.0.and.
     &       esdata(mndata,3).gt.0.0) then
	     sigma(2)=sigma(2)+sdata(mndata,2)/(esdata(mndata,2)**2)
	     esigma(2)=esigma(2)+1.0/(esdata(mndata,2)**2)
	     sigma(3)=sigma(3)+sdata(mndata,3)/(esdata(mndata,3)**2)
	     esigma(3)=esigma(3)+1.0/(esdata(mndata,3)**2)
	  endif
	  if(esdata(mndata,1).le.0.0) then
	    mndata=mndata-1
	    goto 1200
	  else
	    r1(num,mndata)=sdata(mndata,1)
	    re1(num,mndata)=esdata(mndata,1)
	  endif
	  if(sdata(mndata,2).le.0.0.or.sdata(mndata,3).le.0.0) then
	    mndata=mndata-1
	    goto 1200
	  else
	    r2(num,mndata)=sdata(mndata,2)/sdata(mndata,3)
	    re2(num,mndata)=sqrt((sdata(mndata,2)*esdata(mndata,3))
     &	     **2+(sdata(mndata,3)*esdata(mndata,2))**2)
     &	     /(sdata(mndata,3)**2)
	  endif
	  lndata(num)=mndata
 1200	continue
	if(lndata(num).le.0) then
	  write(*,*) 'No data in the period.'
	  write(*,*)
	  num=num-1
	  goto 8900
	endif
*------- Mean value ---------------
	xm1=sigma(1)/esigma(1)
	exm1=1.0/sqrt(esigma(1))
	if(sigma(2).le.0.0.and.sigma(3).le.0.0) then
	   xm2=-99.
	  exm2=-99.
	else
	   xm2=(sigma(2)*esigma(3))/(sigma(3)*esigma(2))
	  exm2=(sqrt((sigma(3)**2)*esigma(2)+(sigma(2)**2)
     &              *esigma(3))*esigma(3))/((sigma(3)**2)*esigma(2))
	endif
	write(*,128) xm1,exm1
  128	format(' Mean value (intensity) =',f8.4,' +-',f8.4)
	write(*,*)
	write(*,129) xm2,exm2
  129	format(' Mean value (color) =',f8.4,' +-',f8.4)
	write(*,*)
*------- Make Comment -----------------------
	write(text(num),600) (ntime0(j),j=1,3),(ntime1(k),k=1,3)
  600	format('  From ',i4,2('/ ',i2),'  To ',i4,2('/ ',i2))
	write(text(num+10),610) lndata(num)
  610	format(8x,'The number of data points = ',i4)
*------- Scale --------------------
	write(*,*) 'Plot all you inputted (y-1/n-0) ?'
*	read(*,*) id
	read(*,'(i2)') id
	write(*,*)
*------- X - jiku -----------------
	xmax1=r1(num,1)
	if(id.eq.1) then
 4001	  do 4000 i41=1,num
 4101	    do 4100 i42=1,lndata(i41)
	      if(xmax1.lt.r1(i41,i42)) xmax1=r1(i41,i42)
 4100	    continue
 4000	  continue
	else
 4201	  do 4200 i43=2,lndata(num)
	    if(xmax1.lt.r1(num,i43)) xmax1=r1(num,i43)
 4200	  continue
	endif
	if(xmax1.le.0.0) xmax1=1.0
*------- Y - jiku -----------------
	ymax1=r2(num,1)
	if(id.eq.1) then
 5001	  do 5000 i51=1,num
 5101	    do 5100 i52=1,lndata(i51)
	      if(ymax1.lt.r2(i51,i52)) ymax1=r2(i51,i52)
 5100	    continue
 5000	  continue
	else
 5201	  do 5200 i53=2,lndata(num)
	    if(ymax1.lt.r2(num,i53)) ymax1=r2(num,i53)
 5200	  continue
	endif
	if(ymax1.le.0.0) ymax1=1.0
*------ Waku wo kaku ------------------------
        if(xmax1.lt.1.0) then
          xmax1=xmax1*100000.
          iflag=1
        endif
        call sutble(xmax1,xmax)
        if(iflag.eq.1) xmax=xmax/100000.
        if(ymax1.lt.1.0) then
          ymax1=ymax1*100000.
          iflag=2
        endif
        call sutble(ymax1,ymax)
        if(iflag.eq.2) ymax=ymax/100000.
	xmin=0.0
	ymin=0.0
	xmax0=xmax
	ymax0=ymax
*------- Plot -------------------------------
 8800	continue
	call icplot
     &  (num,r1,re1,r2,re2,text1,text2,text,cmnt,xmin,xmax,
     &	ymin,ymax,sttmjd,stpmjd,lndata,lday,ldata,ledata,
     &	id,ipage)
*------- Change Scale -------------
	write(*,*)
	write(*,*) 'Change the scales (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 8900
	call scalet(xmin,xmax,ymin,ymax)
	if(xmax.le.0.0) xmax=xmax0
	if(ymax.le.0.0) ymax=ymax0	
	goto 8800
*------- Next -------------------------------
 8900	continue
	write(*,*) 'Plot more (y-1/n-0) ?'
*	read(*,*) iyn
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.eq.1) goto 9000
*------ Return -----------------
  999	continue
	return
	end
******************************************************************
	subroutine icplot
     &  (num,r1,re1,r2,re2,text1,text2,text,cmnt,xmin,xmax,
     &	ymin,ymax,sttmjd,stpmjd,lndata,lday,ldata,ledata,
     &	id,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
        real*4    r1(10,nlim),re1(10,nlim)
	real*4	  r2(10,nlim),re2(10,nlim)
	real*4	  lday(10,200),ldata(10,200),ledata(10,200)
	real*4	  sttmjd(10),stpmjd(10)
	integer	  lndata(10)
 	character cmnt*20,text(30)*50
	character text1*47,text2*47
	integer   id,ipage
*------- Plot start -------------------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
	call pgsch(0.9)
*------- Intensity - Color plot -------------
	call pgvport(0.07,0.6,0.1,0.9)
	call pgscr(25,0.6,0.7,1.0)
	call pgsci(25)
	call pgscf(1)
        call pgwindow(xmin,xmax,ymin,ymax)
        call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------- Plot Coment ------------------------
	call pgsch(0.95)
	call pgscr(17,0.3,1.0,0.1)
	call pgsci(17)
	call pgscf(2)
        call pglabel(text1,text2,'Intensity - Color Diagram')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.83)
	ispec=1
	if(id.eq.1) ispec=num
 6001	do 6000 i6=1,ispec
	  ispp=i6
	  if(id.ne.1) ispp=num
 6101	  do 6100 i61=1,lndata(ispp)
            if(r1(ispp,i61).le.0.0.or.r2(ispp,i61).le.0.0)
     &	      goto 6100
            x0=r1(ispp,i61)
            y0=r2(ispp,i61)
            x1=r1(ispp,i61)-re1(ispp,i61)
            x2=r1(ispp,i61)+re1(ispp,i61)
            y1=r2(ispp,i61)-re2(ispp,i61)
            y2=r2(ispp,i61)+re2(ispp,i61)
	    call pgsci(ispp+6)
	    call pgpoint(1,x0,y0,ispp+1)
            call pgerrx(1,x1,x2,y0,0.0)
            call pgerry(1,x0,y1,y2,0.0)
 6100	  continue
 6000	continue
*------ title ------------------
	call pgscr(18,1.0,0.7,0.7)
	call pgsci(18)
	call pgscf(3)
	call pgsch(0.85)
	x3=xmin+(xmax-xmin)*0.7
	y3=ymin+(ymax-ymin)*0.9
	call pgtext(x3,y3,cmnt(1:20))
*------ Light Curve (1-20 keV) -
	call pgvport(0.68,0.97,0.65,0.9)
	call pgsch(0.78)
*------ X jiku no scale --------
	if(id.eq.1.and.num.gt.1) then
	  xmin2=sttmjd(1)-10.0
	  xmax2=stpmjd(1)+10.0
 6501	  do 6500 i65=2,num
	  if(sttmjd(i65)-10.0.lt.xmin2) xmin2=sttmjd(i65)-10.0
	  if(stpmjd(i65)+10.0.gt.xmax2) xmax2=stpmjd(i65)+10.0
 6500	  continue
	else
	  xmin2=sttmjd(num)-10.0
	  xmax2=stpmjd(num)+10.0
	endif
*------ Y jiku no scale --------
	ymin2=0.0
	ymax2=ldata(num,1)
	if(id.eq.1) then
 8001	  do 8000 i81=1,num
 8101	    do 8100 i82=1,lndata(i81)
	      if(ymax2.lt.ldata(i81,i82)) ymax2=ldata(i81,i82)
 8100	    continue
 8000	  continue
	else
 8201	  do 8200 i83=1,lndata(num)
	    if(ymax2.lt.ldata(num,i83)) ymax2=ldata(num,i83)
 8200	  continue
	endif
*------ Waku wo kaku -----------
	iflag=1
	if(ymax2.le.1.0) then
	  ymax2=ymax2*100000.
	  iflag=4
	endif
	call sutble(ymax2,ymax3)
	if(iflag.eq.4) ymax3=ymax3/100000.
	call pgsci(25)
	call pgscf(1)
	call pgwindow(xmin2,xmax2,ymin2,ymax3)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(19,1.0,0.6,0.2)
	call pgsci(19)
	call pgscf(2)
ce	call pglabel('MJD','Flux (counts/sec/cm**2)',
	call pglabel('MJD','Flux (counts/sec/cm\\u2\\d)',
     &	             'Light Curve (1-20 keV)')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.65)
	ispec=1
	if(id.eq.1) ispec=num
 8301	do 8300 i83=1,ispec
	  ispq=i83
	  if(id.ne.1) ispq=num
 8401	  do 8400 i84=1,lndata(ispq)
	    if(ldata(ispq,i84).le.0.0) goto 8400
	    x6=lday(ispq,i84)
	    y6=ldata(ispq,i84)
	    y7=ldata(ispq,i84)-ledata(ispq,i84)
	    y8=ldata(ispq,i84)+ledata(ispq,i84)
	    call pgsci(ispq+6)
	    call pgpoint(1,x6,y6,ispq+1)
	    call pgerry(1,x6,y7,y8,0.0)
 8400	  continue
 8300	continue
*------ Plot coment ------------
	call pgvport(0.65,1.0,0.0,0.6)
	call pgwindow(0.0,1.0,0.0,1.0)
	call pgsch(0.8)
	call pgscf(1)
	ispr=1
	if(id.eq.1) ispr=num
 7001	do 7000 i7=1,ispr
	  ispt=i7
	  if(id.ne.1) ispt=num
	  call pgsci(ispt+6)
	  x4=0.05
	  y4=1.0-i7*0.12
	  if(y4.lt.0.1) goto 999
	  call pgtext(x4,y4,char(ispt+1)//text(ispt))
	  call pgtext(x4,y4-0.05,text(ispt+10))
 7000	continue
*------- Return -----------------
  999	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
********************************************************
        subroutine dplot8
     & (day,data,edata,ndata,cmnt,xmin,xmax,nbin,ipage)
        parameter (nlim=1600)
        real*4    day(nlim),data(nlim,19),edata(nlim,19)
	real*4	  values(nlim)
ce---	real*4    engch1(16)/1.000,1.674,2.249,2.860,
*    &		             3.479,4.100,4.750,5.338,
*    &		             5.920,7.139,8.314,9.463,
*    &		             10.614,13.034,15.434,17.859/
*	real*4    engch2(16)/1.674,2.249,2.860,3.479,
*    &		             4.100,4.750,5.338,5.920,
*    &		             7.139,8.314,9.463,10.614,
ce---&		             13.034,15.434,17.859,20.284/
	real*4	  sdata(nlim),esdata(nlim)
	integer   nch(2)
	integer   ipage
	character cmnt*20
        character text*45
ce
	real*4	  engch1(16),engch2(16)
        common /engrng/ engch1,engch2
*------- Set Initial -----------
	xmin=0.0
	xmax=0.0
	nbin=10
*------- Input channel ---------
	write(*,*) '*--------- Intensity Distribution ---------*'
	write(*,*)
	ico=0
	iflag=0
*------- distribution ----------
  840	continue
*------- Input A ---------------
	write(*,*) 'Input channel (A1-A2)'
	write(*,*)
	write(*,*) 'A1 ?'
*	read(*,*) nch(1)
	read(*,'(i2)') nch(1)
	write(*,*) 'A2 ?'
*	read(*,*) nch(2)
	read(*,'(i2)') nch(2)
	write(*,*)
*------- Judgement -------------
	if(nch(1).gt.16.or.nch(2).gt.16) goto 840
	if(nch(1).gt.nch(2)) goto 840
*------- plot channel ----------
	write(text,102) engch1(nch(1)),engch2(nch(2))
ce102	format(f5.2,'-',f5.2,' keV counts/sec/cm**2')
  102	format(f5.2,'-',f5.2,' keV counts/sec/cm\\u2\\d')
*------- Clear data ------------
 1101	do 1100 i11=1,nlim
	  sdata(i11)=0.0
	  esdata(i11)=0.0
 1100	continue
*------- Calculate data --------
	ldata=0
 1201	do 1200 i12=1,ndata
	  lch1=nch(1)
	  lch2=nch(2)
	  ldata=ldata+1
	  if(lch1.eq.1.and.lch2.eq.8) then
	    lch1=17
	    lch2=17
	  elseif(lch1.eq.9.and.lch2.eq.16) then
	    lch1=18
	    lch2=18
	  elseif(lch1.eq.1.and.lch2.eq.16) then
	    lch1=19
	    lch2=19
	  endif
	  if(lch1.eq.lch2) then
	    if(edata(i12,lch1).le.0.0) then
	      ldata=ldata-1
	      goto 1200
	    endif
	    sdata(ldata)=data(i12,lch1)
	    esdata(ldata)=edata(i12,lch1)
	  else
 1221	    do 1220 j21=lch1,lch2
	      if(edata(i12,j21).le.0.0) then
		ldata=ldata-1
		goto 1200
	      else
	        sdata(ldata)=sdata(ldata)+data(i12,j21)
	        esdata(ldata)=esdata(ldata)**2
     &	                         +edata(i12,j21)**2
	      endif
 1220	    continue
	    if(esdata(ldata).ge.0.0)
     &	       esdata(ldata)=sqrt(esdata(ldata))
	  endif
 1200	continue
*------- scale -----------------
        if(xmax.eq.0) then
	  xmax=sdata(1)
 3001	  do 3000 i3=1,ldata
	    if(sdata(i3).gt.xmax) xmax=sdata(i3)
 3000	  continue
	endif
	if(nbin.eq.0) nbin=10
*------- set values ------------
 5001	do 5000 i5=1,ldata
	  values(i5)=sdata(i5)
 5000	continue
*------ plot -------------------
	if(ipage.eq.1) call pgbegin(0,'?',1,1)
 6000	continue
        call pgsch(0.9)
	call pgsci(9-ico)
	if(9-ico.eq.7) call pgsci(11)
	call pgscf(1)
	call pghist(ldata,values,xmin,xmax,nbin,iflag)
*------ title ------------------
	call pgsch(0.95)
	call pgsci(7)
	call pgscf(2)
ce	call pglabel('Flux (counts/sec/cm**2)','Occurences'
	call pglabel('Flux (counts/sec/cm\\u2\\d)','Occurences'
     &	 ,'Intensity Distribution of '//cmnt(1:20))
*------- Change Scale -------------
  830	continue
	ico=ico+1
	write(*,*) 'Change the scales (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 820
	write(*,*)
     &   'If you set 0.0, the value will be set automatically.'
	write(*,*) 'X minimum value?'
	read(*,*) xmin
	write(*,*) 'X maximum value?'
	read(*,*) xmax
	write(*,*) 'Bin number?'
	read(*,'(i2)') nbin
	write(*,*)
	goto 6000
*------ next channel -------------
  820	continue
        write(*,*) 'Another one (y-1/n-0)?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 970
	iflag=1
	ipage=0
	goto 840
*------ end -----------------------
  970   continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
        return
        end
********************************************************
	subroutine dplot10
     & (day,data,edata,ndata,cmnt,xmin,xmax,ymin,ymax,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20
	character text*4
	integer   ipage
*------- Graph array set -------
	if(xmax.eq.0.0) xmax=1800.0
	ymin0=ymin
	ymax0=ymax
*------- Plot begin ------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
	call pgsch(0.6)
*------- Iteration start -------
 1001	do 1000 nch=1,16
	ymin=ymin0
	ymax=ymax0
*------- pgvport set -----------
	vx1=(mod(nch-1,4))*0.25
	vy1=(4-int(float(nch-1)/4.0))*0.25
	call pgvport(vx1+0.04,vx1+0.23,vy1-0.2,vy1-0.03)
*------- plot channel ----------
	write(text,'(i2,''ch'')') nch
*------- scale -----------------
	if(ymax.eq.0.0) then
	  ymax1=data(1,nch)
 2001	  do 2000 i2=2,ndata
	    if(ymax1.lt.data(i2,nch)) ymax1=data(i2,nch)
 2000	  continue
	  if(ymax1.le.0.0) ymax1=1.0
	else
	  ymax1=ymax
	endif
*------ Waku wo kaku -----------
	if(ymax.le.0.0) then
	if(ymax1.le.1.0) then
	  ymax1=ymax1*100000.
	  iflag=1
	endif
	call sutble(ymax1,ymax)
	if(iflag.eq.1) ymax=ymax/100000.
	endif
	call pgsci(4)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin,ymax)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgscr(17,1.0,0.5,0.7)
	call pgsci(17)
	call pgscf(2)
	call pglabel('Day','Flux',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.5)
 2201	do 2200 i21=1,ndata
	  if(edata(i21,nch).le.0.0) goto 2200
	  x0=day(i21)
	  y0=data(i21,nch)
	  y1=data(i21,nch)-edata(i21,nch)
	  y2=data(i21,nch)+edata(i21,nch)
	  call pgsci(7)
	  call pgpoint(1,x0,y0,13)
	  call pgerry(1,x0,y1,y2,0.0)
 2200	continue
*------ title ------------------
	call pgsci(3)
	call pgscf(3)
	call pgsch(0.6)
	x1=xmin+(xmax-xmin)*0.78
	y1=ymin+(ymax-ymin)*0.83
	call pgtext(x1,y1,text)
*------ next energy channel ----
 1000	continue
*------ end --------------------
  970	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
************************************************************
	subroutine dplot11
     & (day,data,edata,ndata,cmnt,xmin,xmax,ymin,ymax,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20
	character text*4
	integer   ipage
*------- Graph array set -------
	if(xmax.eq.0.0) xmax=1800.0
	ymin0=ymin
	ymax0=ymax
*------- Plot begin ------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
	call pgsch(0.6)
*------- Iteration start -------
 1001	do 1000 nch=1,16
	ymin=ymin0
	ymax=ymax0
*------- pgvport set -----------
	vx1=(mod(nch-1,4))*0.25
	vy1=(4-int(float(nch-1)/4.0))*0.25
	call pgvport(vx1+0.04,vx1+0.23,vy1-0.2,vy1-0.03)
*------- plot channel ----------
	write(text,'(i2,''ch'')') nch
*------- scale -----------------
	if(ymax.le.0.0) then
	  ymax1=data(1,nch)
 2001	  do 2000 i2=2,ndata
	    if(ymax1.lt.data(i2,nch)) ymax1=data(i2,nch)
 2000	  continue
	  if(ymax1.le.0.0) ymax1=1.0
	else
	  ymax1=ymax
	endif
*------ Kirino ii suuji --------
	iyord=3
	if(ymax.le.0.0) then
	  iloy=int(alog10(ymax1))+1
	  ylog1=float(iloy-1)-alog10(ymax1)
	  if(ylog1.gt.0.0) iloy=iloy-1
	  ymax=iloy
	else
	  ymax=alog10(ymax1)
	endif
	if(ymin.le.0.0) then
	  ymin0=10.**(iloy-iyord)
	  ymin=iloy-iyord
	else
	  ymin0=ymin
	  ymin=alog10(ymin)
	endif
*------ Waku wo kaku -----------
	call pgsci(4)
	call pgscf(1)
	call pgwindow(xmin,xmax,ymin,ymax)
	call pgbox('BCTSN',0.0,0,'BCTSNL',0.0,0)
*------ plot comment -----------
	call pgscr(17,1.0,0.5,0.7)
	call pgsci(17)
	call pgscf(2)
	call pglabel('Day','Flux',' ')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.5)
 2201	do 2200 i21=1,ndata
	  x0=day(i21)
	  if(data(i21,nch).le.0.0) then
	    y0=alog10(ymin0*0.1)
	  else
	    y0=alog10(amax1(data(i21,nch),ymin0))
	  endif
	  y1=alog10(amax1(data(i21,nch)-edata(i21,nch),ymin0))
	  y2=alog10(amax1(data(i21,nch)+edata(i21,nch),ymin0))
	  call pgsci(7)
	  call pgpoint(1,x0,y0,13)
	  call pgerry(1,x0,y1,y2,0.0)
 2200	continue
*------ title ------------------
	call pgsci(3)
	call pgscf(3)
	call pgsch(0.6)
	x1=xmin+(xmax-xmin)*0.78
	y1=ymin+(ymax-ymin)*0.83
	call pgtext(x1,y1,text)
*------ next energy channel ----
 1000	continue
*------ end --------------------
  970	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
****************************************************************
	subroutine spec(day,data,edata,ndata,cmnt,ipage)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	character cmnt*20,text(40)*50
*	real*4    engch(16)/1.744,2.319,2.897,3.472,
*     &		            4.056,4.654,5.222,5.791,
*     &		            6.923,8.041,9.204,10.345,
*     &		            12.718,15.049,17.356,19.869/
ce---	real*4    engch(16)/1.674,2.249,2.860,3.479,
*    &		            4.100,4.750,5.338,5.920,
*    &		            7.139,8.314,9.463,10.614,
ce---&		            13.034,15.434,17.859,20.284/
	dimension ntime0(7),ntime1(7)
	real*8	  mjdhjm,mjdowr,mjday
	real*4    adata(10,16),aedata(10,16)
	real*4	  sdata(10,16),sedata(10,16)
	real*4	  lday(10,nlim),ldata(10,nlim),ledata(10,nlim)
	real*4	  sttmjd(10),stpmjd(10)
	integer   adnum(16)
	integer   isum(10),itch(10)
	integer	  lndata(10)
	integer   ipage
ce
	real*4	  dummy(16),engch(16)
        common /engrng/ dummy,engch
*------- Clear data -------------------
  101	do 100 i11=1,10
  111	  do 110 i12=1,16
	    adata(i11,i12)=0.0
	    aedata(i11,i12)=0.0
	    sdata(i11,i12)=0.0
	    sedata(i11,i12)=0.0
  110	  continue
  100	continue
*------- Program Start ----------------
	write(*,*)
     &	 '*------ Average Spectrum (weighted by error) ------*'
	write(*,*)
	num=0
ce---
 200    continue
        write(*,*) 'Please input data selection mode'
        write(*,*) '1...Day  2...Flux  3...Phase'
        read(*,'(i2)') imode
        write(*,*)
        if (imode.lt.1 .or. imode.gt.3) goto 200
        if (imode.eq.3) then
           write(*,*) 'Input Period (Day)'
           read(*,*) period
           write(*,*)
        endif
        mjdhjm=46796.0d0
        mjdowr=46796.0d0+1800.0d0
ce---
 9000	continue
	num=num+1
	if(num.gt.10) goto 999 
*------- Initialize -------------------
	ymin=0
	ymax=0
	xmin=0
	xmax=0
ce------ Read Day or Discrimination level to add data -------
      if (imode.eq.1) then
 1001   do 1000 i1=4,7
           ntime0(i1)=0
           ntime1(i1)=0
 1000   continue
	write(*,*) 'Input Start day.'
	write(*,*) 
ce	write(*,*) 'Year ?'
ce	read(*,*) ntime0(1)
        write(*,*) 'Year or Day from 1987 Jan 1 ?'
        read(*,*) sday
        if (sday.lt.1987 .or. sday.gt.1992) then
           mjdhjm=dble(sday)+46796.0d0
           write(*,*)
           call mjdate(mjdhjm,ntime0)
           goto 800
        endif
        ntime0(1)=sday
	write(*,*) 'Month ?'
*	read(*,*) ntime0(2)
	read(*,'(i2)') ntime0(2)
	write(*,*) 'Day ?'
*	read(*,*) ntime0(3)
	read(*,'(i2)') ntime0(3)
	write(*,*) 
        call mjulia(ntime0,mjdhjm)
 800    continue
	write(*,*) 'Input Stop day.'
	write(*,*) 
ce	write(*,*) 'Year ?'
ce	read(*,*) ntime1(1)
        write(*,*) 'Year or Day from 1987 Jan 1 ?'
        read(*,*) sday
        if (sday.lt.1987 .or. sday.gt.1992) then
           mjdowr=dble(sday)+46796.0d0
           write(*,*)
           call mjdate(mjdowr,ntime1)
           goto 1300
        endif
        ntime1(1)=sday
	write(*,*) 'Month ?'
*	read(*,*) ntime1(2)
	read(*,'(i2)') ntime1(2)
	write(*,*) 'Day ?'
*	read(*,*) ntime1(3)
	read(*,'(i2)') ntime1(3)
	write(*,*) 
        call mjulia(ntime1,mjdowr)
      elseif (imode.eq.2) then
        write(*,*) 'Input Upper Discrimination.'
        read(*,*) discu
        write(*,*)
        write(*,*) 'Input Lower Discrimination.'
        read(*,*) discl
        write(*,*)
        call mjdate(mjdhjm,ntime0)
        call mjdate(mjdowr,ntime1)
      elseif (imode.eq.3) then
        write(*,*) 'Input Start Phase.'
        read(*,*) stpha
        write(*,*)
        write(*,*) 'Input Stop Phase.'
        read(*,*) sppha
        write(*,*)
        call mjdate(mjdhjm,ntime0)
        call mjdate(mjdowr,ntime1)
      endif
ce---
* 1001	do 1000 i1=4,7
*	  ntime0(i1)=0
*	  ntime1(i1)=0
* 1000	continue
*------- Data no hajime to Owari -------
*	call mjulia(ntime0,mjdhjm)
ce---	call mjulia(ntime1,mjdowr)
 1300 continue
	sttmjd(num)=sngl(mjdhjm)
	stpmjd(num)=sngl(mjdowr)
*------- The number of channels --------
	write(*,*) 'Input the number of the channels ?'
*	read(*,*) itch(num)
	read(*,'(i2)') itch(num)
	write(*,*)
	if(itch(num).le.0.or.itch(num).ge.17) itch(num)=16
*------- Data wo add suru --------------
 1501	do 1500 i15=1,itch(num)
	  adnum(i15)=0
 1500	continue
	lsum=0
	lndata(num)=0
 2001	do 2000 i2=1,ndata
ce	  mjday=dble(day(i2))+46795.5d0
	  mjday=dble(day(i2))+46796.0d0
ce---	  if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 2000
          if(imode.eq.1) then
             if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 2000
          elseif(imode.eq.2) then
             if(data(i2,17).lt.discl .or. data(i2,17).gt.discu)
     &            goto 2000
          elseif(imode.eq.3) then
             phase = day(i2)/period - int(day(i2)/period)
             if(sppha .gt. stpha) then
                if(phase.lt.stpha .or. phase.gt.sppha)
     &               goto 2000
             elseif(sppha .lt. stpha) then
                if(phase.lt.stpha .and. phase.gt.sppha)
     &               goto 2000
             endif
          endif
ce---
	  lsum=lsum+1
*-------- Data ga subete aruka ? -------
 2051	  do 2050 i205=1,itch(num)
	    if(edata(i2,i205).le.0.0) goto 2000
 2050	  continue
*-------- Note! Omomitsuki heikin ------
 2101	  do 2100 i21=1,itch(num)
*	    adata(num,i21)
*     &            =adata(num,i21)+data(i2,i21)
*	    aedata(num,i21)
*     &            =sqrt(aedata(num,i21)**2+edata(i2,i21)**2)
	    adata(num,i21)
     &            =adata(num,i21)+data(i2,i21)/(edata(i2,i21)**2)
	    aedata(num,i21)=aedata(num,i21)+1.0/(edata(i2,i21)**2)
	    adnum(i21)=adnum(i21)+1
 2100	  continue
*-------- Set 1-16 channel data --------
	  lndata(num)=lndata(num)+1
	  lday(num,lndata(num))=sngl(mjday)
	  ldata(num,lndata(num))=data(i2,19)
	  ledata(num,lndata(num))=edata(i2,19)
 2000	continue
*------- No data in the period ---------
	if(lndata(num).eq.0) then
	  write(*,*) 'No Spectrum data in the period.'
	  write(*,*)
	  num=num-1 
	  goto 8000
	endif
*------- Lack of the data --------------
	nlack=lsum-lndata(num)
	write(*,501) nlack
  501	format(1x,'The number of lacked data = ',i4)
	write(*,*)
*------- Data no heikin wo motomeru ----
 3001	do 3000 i3=1,itch(num)
*	  adata(num,i3)=adata(num,i3)/adnum(i3)
*	  aedata(num,i3)=aedata(num,i3)/adnum(i3)
	  adata(num,i3)=adata(num,i3)/aedata(num,i3)
	  aedata(num,i3)=1.0/sqrt(aedata(num,i3))
	  sdata(num,i3)=adata(num,i3)
	  sedata(num,i3)=aedata(num,i3)
 3000	continue
*------- Owarino channel wo matomeru ka? ----
	write(*,*)
     &  'Do you want to synthesize the data for higher channels?'
 3100	continue
	write(*,*) '-------- Input the number --------'
	write(*,*) '        ( No --- input 1 )'
*	read(*,*) isum(num)
	read(*,'(i2)') isum(num)
	write(*,*)
	if(isum(num).le.0.or.isum(num).ge.itch(num)) goto 3100
*------- Synthesize higher channels ---------
	if(isum(num).eq.1) goto 4100
 4001	do 4000 i4=1,isum(num)-1
	  isch=itch(num)-i4
	  adata(num,isch)=adata(num,isch)+adata(num,isch+1)
	  aedata(num,isch)=sqrt(aedata(num,isch)**2+
     &	                   aedata(num,isch+1)**2)
 4000	continue
 4100	continue
	isum(num)=itch(num)-isum(num)+1
*------- Counts/keV ni kanzan ---------------
 3201	do 3200 i32=1,isum(num)
	  if(i32.eq.1) then
	    en1=1.0
	  else
	    en1=engch(i32-1)
	  endif
	  en2=engch(i32)
	  if(i32.eq.isum(num)) en2=engch(itch(num))
	  adata(num,i32)=adata(num,i32)/(en2-en1)
	  aedata(num,i32)=aedata(num,i32)/(en2-en1)
 3200	continue
*------- Make Comment -----------------------
	write(text(num),600) (ntime0(j),j=1,3),(ntime1(k),k=1,3)
  600	format('  From ',i4,2('/ ',i2),'  To ',i4,2('/ ',i2))
	write(text(num+10),610) lndata(num)
  610	format(7x,'The number of added paths = ',i4)
	write(text(num+20),620)(ntime0(j),j=1,3),(ntime1(k),k=1,3)
  620	format('  FROM ',i4,2('/ ',i2),'  TO ',i4,2('/ ',i2))
	write(text(num+30),630) lndata(num)
  630	format(7x,'THE NEMBER OF ADDED PATHS = ',i4)
*------- X jiku no scale --------------------
	xmin=alog10(0.5)
	xmax=alog10(40.0)
*------- Y jiku no scale --------------------
	write(*,*) 'Plot all you inputted (y-1/n-0) ?'
*	read(*,*) id
	read(*,'(i2)') id
	write(*,*)
	ymax1=adata(num,1)
	ymin1=adata(num,1)
	if(id.eq.1) then
 5001	  do 5000 i51=1,num
 5101	    do 5100 i52=1,itch(num)
	      if(ymax1.lt.adata(i51,i52)) ymax1=adata(i51,i52)
	      if(ymin1.gt.adata(i51,i52)) ymin1=adata(i51,i52)
 5100	    continue
 5000	  continue
	else
 5201	  do 5200 i53=1,itch(num)
	    if(ymax1.lt.adata(num,i53)) ymax1=adata(num,i53)
	    if(ymin1.gt.adata(num,i53)) ymin1=adata(num,i53)
 5200	  continue
	endif
	if(ymax1.le.0.0) ymax1=1.0
	if(ymin1.le.0.0) ymin1=0.001
*------- Kirino ii suuji --------------------
	iloy=int(alog10(ymax1))+1
	ylog1=float(iloy-1)-alog10(ymax1)
	if(ylog1.gt.0.0) iloy=iloy-1
	if(ymin.le.0.0) then
	  iloy2=int(alog10(ymin1))-1
	  ymin=iloy2
	  ymin0=10**iloy2
	else
	  ymin0=ymin
	  ymin=int(aint(alog10(ymin)))
	endif
	ymax=iloy
	xmin0=xmin
	xmax0=xmax
	ymin0=ymin
	ymax0=ymax
*------- Graph ni Output suru ---------------
 7000	continue
	call splot
     &  (num,isum,adata,aedata,text,cmnt,xmin,xmax,ymin,ymax,
     &	sttmjd,stpmjd,lndata,lday,ldata,ledata,id,itch,ipage)
*------- Change Scale -------------
	write(*,*)
	write(*,*) 'Change the scales (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.ne.1) goto 8000
	call scalet(xmin,xmax,ymin,ymax)
	if(xmin.le.0.0) then
	  xmin=xmin0
	else
	  xmin=alog10(xmin)
	endif
	if(xmax.le.0.0) then
	  xmax=xmax0
	else
	  xmax=alog10(xmax)
	endif
	if(ymin.le.0.0) then
	  ymin=ymin0
	else
	  ymin=alog10(ymin)
	endif	
	if(ymax.le.0.0) then
	  ymax=ymax0
	else
	  ymax=alog10(ymax)
	endif	
	goto 7000
*------- Next -------------------------------
 8000	continue
	write(*,*) 'Plot more (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.eq.1) goto 9000
*------ Make Spectra file ------------------
	write(*,*)
     &   'Do you want to make a file of all the data inputted'
	write(*,*)
     &   'for spectra fitting (y-1/n-0) ?'
	read(*,'(i2)') iyn
	write(*,*)
	if(iyn.eq.1)
     &    call mksfil(num,sdata,sedata,text,cmnt)
*------ Return -----------------------------
  999	continue
	return
	end
******************************************************************
	subroutine splot
     &  (num,isum,adata,aedata,text,cmnt,xmin,xmax,ymin,ymax,
     &	sttmjd,stpmjd,lndata,lday,ldata,ledata,id,itch,ipage)
	character cmnt*20,text(40)*50
*	real*4    engch(16)/1.744,2.319,2.897,3.472,
*     &		            4.056,4.654,5.222,5.791,
*     &		            6.923,8.041,9.204,10.345,
*     &		            12.718,15.049,17.356,19.869/
ce---	real*4    engch(16)/1.674,2.249,2.860,3.479,
*    &		            4.100,4.750,5.338,5.920,
*    &		            7.139,8.314,9.463,10.614,
ce---&		            13.034,15.434,17.859,20.284/
	real*4    adata(10,16),aedata(10,16)
	real*4	  lday(10,1000),ldata(10,1000),ledata(10,1000)
	real*4	  sttmjd(10),stpmjd(10)
	real*4	  en1,en2
	integer   isum(10),itch(10)
	integer	  lndata(10)
	integer   ipage
ce
	real*4	  dummy(16),engch(16)
        common /engrng/ dummy,engch
*------- Plot start -------------------------
	if(ipage.eq.1) then
	  call pgbegin(0,'?',1,1)
	else
	  call pgpage
	endif
	call pgsch(0.9)
*------- Spectrum plot ----------------------
	call pgvport(0.07,0.6,0.1,0.9)
	call pgscr(25,0.6,0.7,1.0)
	call pgsci(25)
	call pgscf(1)
        call pgwindow(xmin,xmax,ymin,ymax)
        call pgbox('LBCTSN',0.0,0,'LBCTSN',0.0,0)
*------- Plot Coment ------------------------
	call pgsch(0.95)
	call pgsci(3)
	call pgscf(2)
ce	call pglabel('Energy (keV)','Flux (counts/sec/cm**2/keV)'
	call pglabel('Energy (keV)','Flux (counts/sec/cm\\u2\\d/keV)'
     &	             ,'Ginga ASM Spectrum')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.83)
	ymin0=10.0**ymin
	ispec=1
	if(id.eq.1) ispec=num
 6001	do 6000 i6=1,ispec
	  ispp=i6
	  if(id.ne.1) ispp=num
 6101	  do 6100 i61=1,isum(ispp)
	    if(i61.eq.1) then
	      en1=1.0
	    else
	      en1=engch(i61-1)
	    endif
	    en2=engch(i61)
	    if(i61.eq.isum(ispp)) en2=engch(itch(ispp))
	    x0=alog10((en1+en2)*0.5)
	    x1=alog10(en1)
	    x2=alog10(en2)
	    if(adata(ispp,i61).le.0.0) then
	      y0=alog10(ymin0*0.1)
	    else
	      y0=alog10(amax1(adata(ispp,i61),ymin0))
	    endif
	      y1=alog10(amax1(adata(ispp,i61)
     &	                    -aedata(ispp,i61),ymin0))
	      y2=alog10(amax1(adata(ispp,i61)
     &	                    +aedata(ispp,i61),ymin0))
	    call pgsci(ispp+6)
	    call pgpoint(1,x0,y0,ispp+1)
            call pgerrx(1,x1,x2,y0,0.0)
            call pgerry(1,x0,y1,y2,0.0)
 6100	  continue
 6000	continue
*------ title ------------------
	call pgsci(5)
	call pgscf(3)
	call pgsch(0.85)
	x3=xmin+(xmax-xmin)*0.7
	y3=ymin+(ymax-ymin)*0.9
	call pgtext(x3,y3,cmnt(1:20))
*------ Light Curve (1-20 keV) -
	call pgvport(0.68,0.97,0.65,0.9)
	call pgsch(0.78)
*------ X jiku no scale --------
	if(id.eq.1.and.num.gt.1) then
	  xmin2=sttmjd(1)-10.0
	  xmax2=stpmjd(1)+10.0
 6501	  do 6500 i65=2,num
	  if(sttmjd(i65)-10.0.lt.xmin2) xmin2=sttmjd(i65)-10.0
	  if(stpmjd(i65)+10.0.gt.xmax2) xmax2=stpmjd(i65)+10.0
 6500	  continue
	else
	  xmin2=sttmjd(num)-10.0
	  xmax2=stpmjd(num)+10.0
	endif
*------ Y jiku no scale --------
	ymin2=0.0
	ymax2=ldata(num,1)
	if(id.eq.1) then
 8001	  do 8000 i81=1,num
 8101	    do 8100 i82=1,lndata(i81)
	      if(ymax2.lt.ldata(i81,i82)) ymax2=ldata(i81,i82)
 8100	    continue
 8000	  continue
	else
 8201	  do 8200 i83=1,lndata(num)
	    if(ymax2.lt.ldata(num,i83)) ymax2=ldata(num,i83)
 8200	  continue
	endif
*------ Waku wo kaku -----------
	iflag=1
	if(ymax2.le.1.0) then
	  ymax2=ymax2*100000.
	  iflag=4
	endif
	call sutble(ymax2,ymax3)
	if(iflag.eq.4) ymax3=ymax3/100000.
	call pgsci(25)
	call pgscf(1)
	call pgwindow(xmin2,xmax2,ymin2,ymax3)
	call pgbox('BCTSN',0.0,0,'BCTSN',0.0,0)
*------ plot comment -----------
	call pgsci(6)
	call pgscf(2)
ce	call pglabel('MJD','Flux (counts/sec/cm**2)',
	call pglabel('MJD','Flux (counts/sec/cm\\u2\\d)',
     &	             'Light Curve (1-20 keV)')
*------ plot data --------------
	call pgscf(1)
	call pgsch(0.65)
	ispec=1
	if(id.eq.1) ispec=num
 8301	do 8300 i83=1,ispec
	  ispq=i83
	  if(id.ne.1) ispq=num
 8401	  do 8400 i84=1,lndata(ispq)
	    if(ldata(ispq,i84).le.0.0) goto 8400
	    x6=lday(ispq,i84)
	    y6=ldata(ispq,i84)
	    y7=ldata(ispq,i84)-ledata(ispq,i84)
	    y8=ldata(ispq,i84)+ledata(ispq,i84)
	    call pgsci(ispq+6)
	    call pgpoint(1,x6,y6,ispq+1)
	    call pgerry(1,x6,y7,y8,0.0)
 8400	  continue
 8300	continue
*------ Plot coment ------------
	call pgvport(0.65,1.0,0.0,0.6)
	call pgwindow(0.0,1.0,0.0,1.0)
	call pgsch(0.8)
	call pgscf(1)
	ispr=1
	if(id.eq.1) ispr=num
 7001	do 7000 i7=1,ispr
	  ispt=i7
	  if(id.ne.1) ispt=num
	  call pgsci(ispt+6)
	  x4=0.05
	  y4=1.0-i7*0.12
	  if(y4.lt.0.1) goto 999
	  call pgtext(x4,y4,char(ispt+1)//text(ispt))
	  call pgtext(x4,y4-0.05,text(ispt+10))
 7000	continue
*------ next plot --------------
  999	continue
	call pgsch(0.8)
	call pgiden
	write(*,*) 'Change the devices (y-1/n-0) ?'
	read(*,'(i2)') ipage
	if(ipage.eq.1) call pgend
	return
	end
****************************************************************
	subroutine mksfil(num,sdata,sedata,text,cmnt)
	parameter (pp=1.0)
	real*4	  sdata(10,16),sedata(10,16)
ce---	real*4    engch(16)/1.674,2.249,2.860,3.479,
*    &		            4.100,4.750,5.338,5.920,
*    &		            7.139,8.314,9.463,10.614,
ce---&		            13.034,15.434,17.859,20.284/
*	real*4    engch(16)/1.744,2.319,2.897,3.472,
*     &		            4.056,4.654,5.222,5.791,
*     &		            6.923,8.041,9.204,10.345,
*     &		            12.718,15.049,17.356,19.869/
	character cmnt*20,text(40)*50
	character spcfil*30
	logical   exfg1
ce
	real*4	  dummy(16),engch(16)
        common /engrng/ dummy,engch
*-------- Set spec file ----------------
ce	spcfil='spec/'//cmnt(1:20)
	spcfil=cmnt(1:20)
	inquire(file=spcfil,exist=exfg1)
	if(exfg1) then
	  write(*,*)
     &    'The Spectra file is already exists. Please confirm it.'
	  write(*,*)
	  return
	endif
*-------- Open and write in the file ---
	open(11,file=spcfil,status='new')
	write(11,'(a32,127x,a1)')
ce   &                    '  X-ray Spectrum of'//cmnt(1:12),' '
     &                    '  X-ray Spectrum of '//cmnt(1:12),' '
	write(11,'(a23,136x,a1)') '  Before re-arrangement',' '
	write(11,'(a8,141x,a1)') '  Energy',' '
	write(11,'(a8,141x,a1)') '  Counts',' '
	write(11,'(a7,142x,a1)') '  Error',' '
 1001	do 1000 i1=1,num
ce	  write(11,'(''* ADD DATA   16'',2X,A10,2X,a36
ce   &        ,a39,46x,f10.4)') cmnt(2:11),text(i1+20)
	  write(11,'(''* ADD DATA   16'',2X,A12,2X,a36
     &        ,a39,44x,f10.4)') cmnt(1:12),text(i1+20)
     &        ,text(i1+30),pp
          WRITE(11,'(16F10.6)') (engch(j),J=1,16)  
          WRITE(11,'(16F10.4)') (sdata(i1,j),J=1,16)      
          WRITE(11,'(16F10.4)') (sedata(i1,j),J=1,16)
 1000	continue
	close(11)
	write(*,*) 'The file has been completed.'
	write(*,*)
*-------- Return -----------------------
  999	continue
	return
	end                        
****************************************************************
        subroutine mjulia(TIME,MJD)
        INTEGER*4  TIME(7)
        REAL*8     MJD
C
C  DATE AND TIME ----> MJD
C
C     INPUT DATA; TIME(7)
C                  YEAR, MONTH, DAY, HOUR, MINUTE, SEC, MILISEC
C                  IF YEAR < 100  , YEAR IS ASSUME TO BE 1900'S
C     OUT PUT DATA; MJD
C                  MJD=JD-2400000.5
C            DAY OF WEEK = MOD(MJD+3,7)  (0=SUN, 1=MON, ETC)
C
        REAL*8 Y,M,D
        IF(TIME(2).GT.2) THEN
          Y=DFLOAT(TIME(1))
          M=DFLOAT(TIME(2))
         ELSE
          Y=DFLOAT(TIME(1)-1)
          M=DFLOAT(TIME(2)+12)
        ENDIF
        IF(TIME(1).LT.100) Y=Y+1900.0
        D=DFLOAT(TIME(3))+
     &    (DFLOAT(TIME(4))+(DFLOAT(TIME(5))+
     &    (DFLOAT(TIME(6))+DFLOAT(TIME(7))/1000.0)/60.0 )/60.0) 
     &	   /24.0
        MJD=DINT(Y*365.25)+DINT(Y/400.)-DINT(Y/100.)
     &       +DINT((M-2.)*30.59) + D -678912.0D0
        RETURN
        END
****************************************************************
        subroutine mjdate(mjd,time)
        real*8    mjd
        INTEGER   TIME(7)
C
C  MJD ---> DATE AND TIME
C                           87/5/18 MITSUDASAN NO SHIJIDE HENKOU
C     INPUT DATA; MJD                                BY T.DOTANI
C                  MJD=JD-2400000.5
C     OUTPUT DATA; TIME(7)
C                  YEAR, MONTH, DAY, HOUR, MINUTE, SEC, MILISEC
C
        REAL*8 MJD0,DAY,MJD1,MJD2
        INTEGER*4 MONTH(12)/31,28,31,30,31,30,31,31,30,31,30,31/
        INTEGER   J,TIME2(7)/0,1,1,0,0,0,0/
        integer   flag

        MJD1=MJD+5.7870370D-9
        TIME(1)=INT(MJD1*0.0027379093D0+1858.877D0)
        TIME(2)=1
        TIME(3)=1
        DO 10 J=4,7
          TIME(J)=0
   10   CONTINUE
        FLAG=0
   30   continue
          TIME2(1)=TIME(1)+1
          CALL MJULIA(TIME,MJD0)
          CALL MJULIA(TIME2,MJD2)
          IF(MJD1.LT.MJD0) THEN
            TIME(1)=TIME(1)-1
          ELSEIF(MJD1.GE.MJD2) THEN
            TIME(1)=TIME(1)+1
          ELSE
            FLAG=1
            goto 32
          ENDIF
          goto 30
   32   CONTINUE
        IF(MOD(TIME(1),4).EQ.0.AND.MOD(TIME(1),100).NE.0) THEN
          MONTH(2)=29
         ELSE
          MONTH(2)=28
        ENDIF
        DAY=MJD1-MJD0+1.0D0
        J=1
      	DO 20 WHILE(INT(DAY).GT.MONTH(J))
       	 DAY=DAY-DFLOAT(MONTH(J))
         J=J+1
20    	CONTINUE
      	TIME(2)=J
      	IF(DAY.LT.33.0) THEN
       	 TIME(3)=DAY
         DAY=(DAY-TIME(3))*24
         TIME(4)=DAY
         DAY=(DAY-TIME(4))*60
         TIME(5)=DAY
         DAY=(DAY-TIME(5))*60
         TIME(6)=DAY
         TIME(7)=(DAY-TIME(6))*1000
     	 ENDIF
        RETURN
        END
***************************************************************
	subroutine wrtfil(fname,day,data,edata,ndata,cmnt)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	real*4	  sday(nlim),sdata(nlim,19),esdata(nlim,19)
	real*8	  mjdhjm,mjdowr,mjday
	integer	  ntime0(7),ntime1(7)
	character cmnt*20
	character fname*10
	character fname1*20,fname2*20,fname3*20
	logical   exfg1,exfg2,exfg3
*------- Clear data --------------
 6001	do 6000 i6=1,nlim
	  sday(i6)=0.0
 6101	  do 6100 i61=1,19
	    sdata(i6,i61)=0.0
	    esdata(i6,i61)=0.0
 6100	  continue
 6000	continue
*------- file name ---------------
	len=10
 1001	do 1000 i11=1,10
	  if(fname(i11:i11).eq.' ') then
	    len=i11-1
	    goto 1010
	  endif
 1000	continue
 1010	continue
ce---	write(fname1,'(''macfile/'',a,''.a'')') fname(1:len)
*	write(fname2,'(''macfile/'',a,''.b'')') fname(1:len)
ce---	write(fname3,'(''macfile/'',a,''.c'')') fname(1:len)
	write(fname1,'(a,''.a'')') fname(1:len)
	write(fname2,'(a,''.b'')') fname(1:len)
	write(fname3,'(a,''.c'')') fname(1:len)
*------- Search file -------------
	inquire(file=fname1,exist=exfg1)
	if(exfg1)then
	  write(*,120) fname1
	  write(*,*)
	  return
	endif
	inquire(file=fname2,exist=exfg2)
	if(exfg2)then
	  write(*,120) fname2
	  write(*,*)
	  return
	endif
	inquire(file=fname3,exist=exfg3)
	if(exfg3)then
	  write(*,120) fname3
	  write(*,*)
	  return
	endif
  120	format('The file ',a20,
     &         'already exists. Please confirm it.')
*------- hajime to owari no hi ---
	ntime0(1)=87
	ntime0(2)=1
	ntime0(3)=1
	write(*,*) 'Input Stop day.'
	write(*,*) 
	write(*,*) 'Year ?'
*	read(*,*) ntime1(1)
	read(*,'(i2)') ntime1(1)
	write(*,*) 'Month ?'
*	read(*,*) ntime1(2)
	read(*,'(i2)') ntime1(2)
	write(*,*) 'Day ?'
*	read(*,*) ntime1(3)
	read(*,'(i2)') ntime1(3)
	write(*,*) 
 7001	do 7000 i7=4,7
	  ntime0(i7)=0
	  ntime1(i7)=0
 7000	continue
	call mjulia(ntime0,mjdhjm)
	call mjulia(ntime1,mjdowr)
*------- Make data ----------------------------
	ldata=0
 6201	do 6200 i62=1,ndata
ce	  mjday=dble(day(i62))+46795.5d0
	  mjday=dble(day(i62))+46796.0d0
	  if(mjday.lt.mjdhjm.or.mjday.gt.mjdowr) goto 6200
	  ldata=ldata+1
	  sday(ldata)=day(i62)
	  icount=0
 6301	  do 6300 i63=1,19
	    if(edata(i62,i63).le.0.0) then
	      sdata(ldata,i63)=-9.9
	      esdata(ldata,i63)=-9.9
	      icount=icount+1
	    else
	      sdata(ldata,i63)=data(i62,i63)
	      esdata(ldata,i63)=edata(i62,i63)
	    endif
 6300	  continue
	  if(icount.ge.19) ldata=ldata-1
 6200	continue
*------- Output data for 2 channels -----------
	open(66,file=fname1,status='new')
	write(66,'(a)') ' ***** Result of ASM scanning *****'
	write(66,'('' Target = '',a10)') cmnt(1:20)
	write(66,100) (ntime0(i),i=1,3),(ntime1(j),j=1,3)
  100	format(' From ',i4,2('/ ',i2),'  To ',i4,2('/ ',i2))
	write(66,*) ' *------------- 1-16 ch data -------------*'
	write(66,*) 'Num      Day    1-6keV   Error  6-20keV',
     &		    '  Error  1-20keV  Error'
 5001	do 5000 i5=1,ldata
	  write(66,300)
     &	    i5,sday(i5),(sdata(i5,m),esdata(i5,m),m=17,19)
  300	  format(i4,f11.3,6(1x,f7.3))
 5000	continue
	close(66)
*------- Output data for 1-8 channels ----------
	open(66,file=fname2,status='new')
	write(66,'(a)') ' ***** Result of ASM scanning *****'
	write(66,'('' Target = '',a10)') cmnt(1:20)
	write(66,100) (ntime0(i),i=1,3),(ntime1(j),j=1,3)
	write(66,*) ' *------------- 1-8 ch data -------------*'
	write(66,310) (mn,mn=1,8)
  310	format('Num Day',8(1x,i2,'ch Error')) 
 8001	do 8000 i8=1,ldata
	  write(66,320)
     &	    i8,sday(i8),(sdata(i8,m),esdata(i8,m),m=1,8)
  320	  format(i4,f11.3,16(1x,f6.3))
 8000	continue
	close(66)
*------- Output data for 9-16 channels ---------
	open(66,file=fname3,status='new')
	write(66,'(a)') ' ***** Result of ASM scanning *****'
	write(66,'('' Target = '',a10)') cmnt(1:20)
	write(66,100) (ntime0(i),i=1,3),(ntime1(j),j=1,3)
	write(66,*) ' *------------- 9-16 ch data -------------*'
	write(66,310) (mn,mn=9,16)
 8101	do 8100 i81=1,ldata
	  write(66,320)
     &	    i81,sday(i81),(sdata(i81,m),esdata(i81,m),m=9,16)
 8100	continue
	close(66)
*------- Return-------------------
  999	continue
	write(*,*) 'The file has been completed.'
	write(*,*)
	return
	end
***************************************************************
	subroutine matome(day,data,edata,ndata)
	parameter (nlim=1600)
	real*4    day(nlim),data(nlim,19),edata(nlim,19)
	real*4	  bday(nlim),bdata(nlim,19),ebdata(nlim,19)
	real*4	  gdata(19),egdata(19)
	real*4	  maxd
*------- Clear data -----------
 1001	do 1000 i1=1,ndata
 1101	  do 1100 i11=1,19
	    bdata(i1,i11)=0.0
	    ebdata(i1,i11)=0.0
 1100	  continue
 1000	continue
*------- Input maotme day -----
 1200	continue
	write(*,*) 'Input duration accumulated the data.'
	write(*,*)
*	read(*,*) nday
	read(*,'(i2)') nday
	write(*,*)
*------- Judge ----------------
	if(nday.gt.1800) goto 1200
*------- max data -------------
	maxd=day(ndata)
 2001	do 2000 i2=1,ndata-1
	  if(maxd.lt.day(i2)) maxd=day(i2)
 2000	continue
	nn=maxd/float(nday)+3
*------- accumulate -----------
	mdata=0
 3001	do 3000 i3=1,nn
	  stt=float(nday*(i3-1))
	  stp=float(nday*i3)
	  gday=0.0
 3401	  do 3400 i34=1,19
	    gdata(i34)=0.0
	    egdata(i34)=0.0
 3400	  continue
	  mdata=mdata+1
	  ns=0
 3101	  do 3100 i31=1,ndata
	    if(day(i31).lt.stt.or.day(i31).gt.stp) goto 3100
	    gday=gday+day(i31)
	    ns=ns+1
 3201	    do 3200 i32=1,19
	      if(data(i31,i32).le.0.0) goto 3200
	      gdata(i32)=gdata(i32)+data(i31,i32)
     &		                 /(edata(i31,i32)**2)
	      egdata(i32)=egdata(i32)+1.0/(edata(i31,i32)**2)
 3200	    continue
 3100	  continue
	  if(ns.le.0) then
	    mdata=mdata-1
	    goto 3000
	  else
	    bday(mdata)=gday/float(ns)
 3301	    do 3300 i33=1,19
	      if(gdata(i33).le.0.0) then
		bdata(mdata,i33)=-99.
		ebdata(mdata,i33)=-99.
	      else
	        bdata(mdata,i33)=gdata(i33)/egdata(i33)
	        ebdata(mdata,i33)=sqrt(1.0/egdata(i33))
	      endif
 3300	    continue
	  endif
 3000	continue
*------- interchange data -------
 4001	do 4000 i41=1,ndata
	  day(i41)=bday(i41)
 4101	  do 4100 i42=1,19
	    data(i41,i42)=bdata(i41,i42)
	    edata(i41,i42)=ebdata(i41,i42)
 4100	  continue
 4000	continue
*------- return -----------------
  999	continue
	return
	end
