      program convf2t
      use dflib
      use dfport
      implicit real*8 (a-h,o-z)
      parameter (md=9999)
      integer status,unit,readwrite,blocksize,naxes(2),nfound
      character*60 fname,fnout,fnouttemp
      character*80 comment
      character*10 odate,ut,ra,dec,object
      real*4 sltwid,exptime
      logical LOREFR
      character*60 fnrefspec,labref
      common/reflab/LOREFR,fnrefspec,labref
      common/wrange/wmne,wmxe,jbl,jbr,lolimit
      logical lolimit
      dimension wlasel(md),specsel(md)
c
      common/orderinfo/numord,nump(99),wlmul(md,99),spmul(md,99),
     & nofoc,nobest,novtot,novis(99),xfrac(99)
      common/dummyorder/npx,nmar
c
      data itemp/0/
      dimension np(99)
      real*4 x(md,99),y(md,99),xl(2),yl(2)
      real*4 wcen,wid,wmin,wmax,smin,smax,shift
      character*133 label
      character*120 fnselist
      character*80 line
      character*3 chnum,czu
      integer*2 narg1,narg2
      character*80 buf
      logical back
c----------------------------
      lorefr=.false.
      npx=9999
      nmar=0
c----------------------------

      narg1=1
      call getarg(narg1,buf)
      read(buf,'(a)') fname
c      write(*,'(a)') 'arg1=',fname

      fnouttemp=''
      iffo=0
      narg2=2
      call getarg(narg2,buf)
      read(buf,'(a)',err=19,end=19) fnouttemp
c      write(*,'(a)') 'arg2=',fnouttemp
      if(fnouttemp.eq.fname) go to 19
      if(fnouttemp.eq.'') go to 19
      if(index(fnouttemp,'.txt').le.0) then
        write(*,*) 'why not extention "*.txt"?-->'//trim(fnouttemp)
        stop
      end if
      iffo=1
   19 continue

      ift=index(fname,'.fit',back=.true.)
      if(ift.le.0) then
      write(*,*) 
     & 'Convf2t: Input filename should have extention "*.fit(s)"'
      stop
      end if

c      write(*,'(a)') 'Specified FITS file to convert: '//trim(fname)
C     open the FITS file
c      readwrite=0
c 3    call ftopen(98,fname,readwrite,blocksize,status)
c      call ftgkys(98,'OBJECT  ',object,comment,status)
c      call ftgkys(98,'DATE-OBS',odate,comment,status)
c      call ftgkys(98,'UT      ',ut,comment,status)
c      call ftgkys(98,'RA      ',ra,comment,status)
c      call ftgkys(98,'DEC     ',dec,comment,status)
c      call ftgkye(98,'SLT-WID ',sltwid,comment,status)
c      call ftgkye(98,'EXPTIME ',exptime,comment,status)
cC     close the file and free the unit number
c 7    call ftclos(98, status)
c      call ftfiou(98, status)
      call checkfit(fname,ifbl,norder)
      call fitsread(98,fname,ifbl,norder,nptsel,wlasel,specsel)
c      write(21,100) fname(ift-16:ift+4),object,ra,dec,odate,ut,
c     % sltwid,exptime,numord
c  100 format(a21,5(2x,a10),f9.3,f7.0,2x,i5)
c      do j=1,numord
c      write(21,'(i4,2(f11.3,f9.5))') 
c     & nump(j),wlmul(1,j),spmul(1,j),wlmul(nump(j),j),spmul(nump(j),j)
c      end do

      if(iffo.eq.1) then
      fnout=fnouttemp
      go to 88
      end if

c     now setting default fnout
      if(numord.le.1) then
      czu='_ss'
      else
      czu='_ms'
      end if
      ify=index(fname,'\',back=.true.)
      if(ify.gt.0) then
      fnout=fname(ify+1:ift-1)//czu//'.txt'
      else
      fnout=fname(1:ift-1)//czu//'.txt'
      end if

   88 continue
c      write(*,'(a)') 
c     & 'Now writing the data into text file --> '//trim(fnout)
      open (unit=97,file=fnout)
      write(97,'(a)') 
     & '% Converted from fits-formatted file:'//trim(fname)
      write(97,'(''% number of orders:'',i4)') numord
      write(97,'(a)') 
     & '% Structure of data points and wavelength regions:'  
      do j=1,numord
        if(wlmul(1,j).gt.wlmul(nump(j),j)) then
          write(97,'(''% '',i3,i6,2f11.3)') 
     &    j, nump(j),wlmul(nump(j),j),wlmul(1,j)
        else
          write(97,'(''% '',i3,i6,2f11.3)') 
     &    j, nump(j),wlmul(1,j),wlmul(nump(j),j)
        end if
      end do
      write(97,*)
      do j=1,numord
        if(wlmul(1,j).gt.wlmul(nump(j),j)) then
          do n=nump(j),1,-1
            write(97,'(f15.6,e20.10)') wlmul(n,j),spmul(n,j) 
          end do
        else
          do n=1,nump(j)
            write(97,'(f15.6,e20.10)') wlmul(n,j),spmul(n,j) 
          end do
        end if
        if(j.lt.numord) write(97,'(f15.6,e20.10)') 0., 0.
      end do
      close (97)
      stop
      end

      subroutine checkfit(fnobspec,ifbl,norder)
      implicit real*8 (a-h,o-z)
      character*60 fnobspec
      common/wrange/wmne,wmxe,jbl,jbr,lolimit
      logical lolimit
c
      logical LOREFR
      character*60 fnrefspec,labref
      common/reflab/LOREFR,fnrefspec,labref
c
c     norder is returned
c     norder=-999  this is not a fits-formatted file (treated as text file)
c     norder=0     single-spectrum fits-formatted file
c     norder>0     multi-spectra fits-formatted file --> 'norder' selected
c     ifbl is returned
c     |ifbl| indicates the column number of '[' character
c     ifbl>0 ---> [order,wavelength-points]
c     ifbl<0 ---> [wavelength-points,order]
c

c      ifpfit =scan(fnobspec,'.fit')
      call lablen(fnobspec,60,ns,ne)
      do i=ns,ne-3
      if(fnobspec(i:i+3).eq.'.fit'.or.fnobspec(i:i+3).eq.'.FIT') 
     & go to 76
      end do
      ifpfit=0
      go to 77
   76 ifpfit=i
   77 continue

      ifbl=scan(fnobspec,'[')
      ifbr=scan(fnobspec,']')
      ifcomma=scan(fnobspec,',')
      ifstar=scan(fnobspec,'*')
c#
      IF(.NOT.LOREFR) THEN
      wmne=0.
      wmxe=1.e30
      lolimit=.false.
      if(ifbl*ifbr.ne.0.and.ifstar.eq.0) then
      lolimit=.true.
      jbl=ifbl
      jbr=ifbr
      read(fnobspec(jbl+1:jbr-1),*) wmne,wmxe
      end if
      END IF
c#
c
      if(ifpfit.eq.0) then
      norder=-999
      ifbl=0
      return
      end if

      if(ifbl*ifbr*ifcomma*ifstar.eq.0) then
      norder=0
      return
      end if
c
      if(ifcomma.lt.ifstar) then
      read(fnobspec(ifbl+1:ifcomma-1),*) norder
      ifbl=+iabs(ifbl)
      else
      read(fnobspec(ifcomma+1:ifbr-1),*) norder
      ifbl=-iabs(ifbl)
      end if
      return
      end


      subroutine fitsread(nn,fnamet,ifbl,norder,nptsel,wlasel,specsel)
      use dflib
      implicit real*8 (a-h,o-z)
      parameter (md=9999)
c#
      common/wrange/wmne,wmxe,jbl,jbr,lolimit
      logical lolimit
      dimension wlasel(md),specsel(md)
      real*8,dimension(:),allocatable :: wla,spec
c      real*4,dimension(:),allocatable :: buffer
c      real*8 wla(md),spec(md)
      real*4 buffer(md*100)
      character*60 fnamet,fname
c#
      integer status,unit,readwrite,blocksize,naxes(2),nfound
      integer group,firstpix,nbuffer,npixels,i
      integer npoint,npt,norder
      real*4 nullval
      logical anynull,loexist
      character filename*80
      character comment*80,ctype1*8,waxmap01*80
      real*4 crval1,cdelt1,crpix1,cd1_1
      real*4 ltv1, ltm1_1

      integer jdcflag,nwcsdim
      character word8*8,fmt*80,waxmap*80,segment(999)*80,string*2048
      dimension ns(999),ne(999)
      character linetemp*80,card(9999)*80
      dimension coef(99),x(99), npln(99),noff(99)
c
      logical LOREFR
      character*60 fnrefspec,labref
      common/reflab/LOREFR,fnrefspec,labref
c
c
      common/orderinfo/numord,nump(99),wlmul(md,99),spmul(md,99),
     & nofoc,nobest,novtot,novis(99),xfrac(99)
      common/dummyorder/npx,nmar
c
      exp10(x)=exp(x*2.302585092994e0)
c
c     checking the file name
      IF(.NOT.LOREFR) numord=1
      if(ifbl.eq.0) then
      fname=fnamet
      else
      fname=fnamet(1:iabs(ifbl)-1)
      end if
      inquire(file=fname,exist=loexist)
      if(.not.loexist) then
      write(*,*) 'specified spectrum file does not exist!'
      stop
      end if

c     open the FITS file
      unit=nn
      readwrite=0
   3  call ftopen(unit,fname,readwrite,blocksize,status)

c     listing-up header data
      kk=0
      do key=1,9999
        call ftgrec(unit,key,linetemp,status)
        call lablen(linetemp,80,nns,nne)
        if(linetemp(1:3).eq.'END') go to 79
        if (nns.gt.nne) cycle
        kk=kk+1
        card(kk)=linetemp
c      write(*,'(i4,2x,a)') kk,trim(card(kk))
      end do
   79 nline=kk

C     determine the size of the image and the format type

      inaxis=0
      inaxis1=0
      inaxix2=0
      ictype1=0
      do k=1,nline
        if(index(card(k),'NAXIS   ').gt.0) inaxis=k
        if(index(card(k),'NAXIS1  ').gt.0) inaxis1=k
        if(index(card(k),'NAXIS2  ').gt.0) inaxis2=k
        if(index(card(k),'CTYPE1  ').gt.0) ictype1=k
      end do
      read(card(inaxis)(10:),*) ndim
      read(card(inaxis1)(10:),*) naxes(1)
      if(ndim.ge.2) then
        read(card(inaxis2)(10:),*) naxes(2)
      else
        naxes(2)=0
      end if
      read(card(ictype1)(10:),*) ctype1

c      write(*,*) 'ndim    =',ndim
c      write(*,*) 'naxes(1)=',naxes(1)
c      write(*,*) 'naxes(2)=',naxes(2)
c      write(*,'(''CTYPE1:'',a8)') ctype1

      if(naxes(1).gt.md*99) then
        write(*,*) 'naxes(1) > md*99 !'
        stop
      end if
      if(naxes(2).gt.99) then
        write(*,*) 'naxes(1) > 99 !'
        stop
      end if


c     check dispersion-related parameters

      iwcsdim=0
      iltm1_1=0
      iltv1=0
      iwaxmap01=0
      do k=1,nline
        if(index(card(k),'WCSDIM  ').gt.0) iwcsdim=k
        if(index(card(k),'LTM1_1  ').gt.0) iltm1_1=k
        if(index(card(k),'LTV1    ').gt.0) iltv1=k
        if(index(card(k),'WAXMAP01').gt.0) iwaxmap01=k
      end do

      if(iwcsdim.gt.0) then
        read(card(iwcsdim)(10:),*) nwcsdim
      else
        nwcsdim=ndim
      end if

      if(iltv1.gt.0) then
        read(card(iltv1)(10:),*) ltv1
      else
        ltv1=0.
      end if

      if(iltm1_1.gt.0) then
        read(card(iltm1_1)(10:),*) ltm1_1
      else
        ltm1_1=1.
      end if

      if(iwaxmap01.gt.0) then
        read(card(iwaxmap01)(10:),*) waxmap01
      else
        waxmap01=''
      end if

c      write(*,'(''WCSDIM  :'',i4)') nwcsdim
c      write(*,'(''LTV1    :'',f8.5)') ltv1
c      write(*,'(''LTM1_1  :'',f8.5)') ltm1_1
c      write(*,'(''WAXMAP01:'',(a))') trim(waxmap01)


cc


      if(ctype1.eq.'MULTISPE') go to 1001

c     non-MULTISPEC case (e.g.,LINEAR dispersion)

      icrval1=0
      icrpix1=0
      icdelt1=0
      icd1_1=0
      idcflag=0
      do k=1,nline
        if(index(card(k),'CRVAL1').gt.0) icrval1=k
        if(index(card(k),'CRPIX1').gt.0) icrpix1=k
        if(index(card(k),'CDELT1').gt.0) icdelt1=k
        if(index(card(k),'CD1_1').gt.0) icd1_1=k
        if(index(card(k),'DC-FLAG').gt.0) idcflag=k
      end do

      if(icrval1*icrpix1.eq.0) then
        write(*,*) 'Either CRVAL1 or CRPIX1 not found'
        write(*,*) 'despite the non-MULTISPEC (single spectrum) mode...'
        write(*,*) 'stopped!'
      stop
      end if

      if(icdelt1.eq.0.and.icd1_1.eq.0) then
        write(*,*) 'Both CDELT1 and CD1_1 not found'
        write(*,*) 'despite the non-MULTISPEC (single spectrum) mode...'
        write(*,*) 'stopped!'
      stop
      end if      

      read(card(icrval1)(10:80),*) crval1
      read(card(icrpix1)(10:80),*) crpix1
      if(icdelt1.ne.0) read(card(icdelt1)(10:80),*) cdelt1
      if(icd1_1.ne.0)  read(card(icd1_1)(10:80),*) cd1_1
      if(icd1_1.eq.0) cd1_1=cdelt1



      if(idcflag.gt.0) then
        read(card(idcflag)(10:),*) jdcflag
      else
        jdcflag=0
      end if

c      write(*,'(''CRVAL1:'',f8.2)') crval1
c      write(*,'(''CRPIX1:'',f8.2)') crpix1
c      write(*,'(''CDELT1:'',f8.5)') cdelt1
c      write(*,'(''CD1_1 :'',f8.5)') cd1_1
c      write(*,'(''DC-FLAG  :'',i5)') jdcflag

      npoint=naxes(1)
      allocate(wla(npoint))
      allocate(spec(npoint))
      
      do l=1,npoint
        wla(l)=crval1+cd1_1*(l-crpix1)
        if(jdcflag.eq.1) wla(l)=exp10(wla(l))
      end do
      npt=npoint
      go to 1003

 1001 continue
ccc   MULTISPE
c      if(norder.eq.0) then
c        call ftgkys(unit,'WAXMAP01',waxmap,comment,status)
c        read(waxmap,*) idum,idum,idum,iposm
c        ipos=iposm+1
c      else
c        ipos=norder
c      end if

      ifdash=0
      ifinv=0

      do 450 i=1,999
        if(i.le.9) then
        fmt='(''WAT2_00'',i1)'
        else if (i.le.99) then
        fmt='(''WAT2_0'',i2)'
        else if (i.le.999) then
        fmt='(''WAT2_'',i3)'
        end if
        write(word8,fmt) i
        do 451 k=1,nline
        if(index(card(k),word8).gt.0) go to 456
  451   continue
        go to 459
  456   segment(i)=card(k)(12:79)
        call lablen(segment(i),80,ns(i),ne(i))
        if(segment(i)(1:1).eq.' ') ns(i)=ns(i)-1
        if(segment(i)(68:68).eq.' ') ne(i)=ne(i)+1
  450 continue
  459 numstr=i-1

      do 8000 j=1,99
      ipos=j
      ic=0
      do i=1,numstr
        do 43 n=ns(i),ne(i)
        if(segment(i)(n:n).ne.'"') go to 43
        ic=ic+1
        if(ic.eq.2*ipos-1) then
          ii1=i
          n1=n
        end if
        if(ic.eq.2*ipos) then
          ii2=i
          n2=n
        go to 44
        end if
   43   continue
      end do
   44 continue
      iii=ii2-ii1+1
      if(iii.gt.99) then
      write(*,*)'Too long WAT2 label at an order (> 99 lines!)-->', iii
      stop
      end if
      if(iii.eq.1) then
        string=segment(ii1)(n1:n2)
      else if (iii.eq.2) then
        string=segment(ii1)(n1:ne(ii1))//segment(ii2)(ns(ii2):n2)
      else if (iii.ge.3) then
        string=segment(ii1)(n1:ne(ii1))
        nnn=ne(ii1)-n1+1
        do iix=ii1+1,ii2-1
          string=string(1:nnn)//segment(iix)(ns(iix):ne(iix))
          nnn=nnn+ne(iix)-ns(iix)+1
        end do
        string=string(1:nnn)//segment(ii2)(ns(ii2):n2)
      end if
      call lablen(string,2048,nt1,nt2)
c      write(*,*) 'nt1,nt2 ',nt1,nt2
c      write(*,'(a)') '|'//string(nt1:nt2)//'|'
      read(string(nt1+1:nt2-1),*) iord,iduma,idumb,crval1,cdelt1,
     & npdum,voc,apb1,apb2
c
      iordp=iord
c#     special treatment for HDS header with dash
      if(j.eq.1.and.iord.eq.iduma) then
        iordabs1=iord
        ifdash=1
      end if
      if(ifdash.eq.1)iordp=iord-iordabs1+1
c#
c%    inverse order case (e.g., old Xinglong CES spectrum)
      if(j.eq.1.and.iord.ne.1.and.iord.eq.naxes(2)) then
        ifinv=1
        iordmax=iord
      end if
      if(ifinv.eq.1) iordp=iordmax-iord+1
c%
c----------------------------------------------------------------------
      if(idumb.eq.0.or.idumb.eq.1) then
        go to 801
      else if (idumb.eq.2) then
        go to 802
      else
        write(*,*) 'Why not idumb=0 or 1 or 2 ?! idumb--> ',idumb
        stop
      end if
  801 continue
c     idumb=0 (for the linear case) or idumb=1 (for the log-linear case)

c      write(*,'(''CRVAL1:'',f8.2)') crval1
c      write(*,'(''CDELT1:'',f8.5)') cdelt1
      w0=crval1
      dw=cdelt1
c      if(ifbl.gt.0) then
c      write(*,*) 'ifbl should be <0'
c      write(*,*) 'transpose row/column as [wavelength-points,order]'
c      stop
c      end if
c      IF(.NOT.LOREFR) THEN
      nump(j)=naxes(1)
      if(idumb.eq.0) then
c     linear (idumb=0)
        do k=1,nump(j)
          wlmul(k,j)=(w0+(k-1)*dw)/(1.+voc)
        end do
      else
c     log-linear (idumb=1)
        do k=1,nump(j)
          wlmul(k,j)=exp10((w0+(k-1)*dw))/(1.+voc)
        end do
      end if
c      write(*,*) j,nump(j),wlmul(1,j),wlmul(nump(j),j)
c      END IF
      go to 7777
  802 continue
c     idumb=2 (for the non-linear dispersion case )
        read(string(nt1+1:nt2-1),*) iord,iduma,idumb,crval1,cdelt1,
     &  npdum,voc,apb1,apb2,wt,offset,iftype
      
      nump(j)=naxes(1)
c
      if(iftype.eq.1) then
        go to 701
      else if (iftype.eq.2) then
        go to 702
      else if (iftype.eq.3) then
        go to 703
      else if (iftype.eq.4) then
        go to 704
cc      else if (iftype.eq.5) then
cc        go to 705
cc      else if (iftype.eq.6) then
cc        go to 706
      else
        write(*,*) 
     &  'convf2t of this version does not support iftype -->',iftype
        stop
      end if

c==== chebyshev case
  701 continue
      read(string(nt1+1:nt2-1),*) iord,iduma,idumb,crval1,cdelt1,
     &  npdum,voc,apb1,apb2,wt,offset,iftype,nodr,pmin,pmax,
     & (coef(l),l=1,nodr)
      do 171 k=1,nump(j)
      p=(k-ltv1)/ltm1_1
      en=(p-(pmax+pmin)/2.)/((pmax-pmin)/2.)
      x(1)=1.
      x(2)=en
      do l=3,nodr
      x(l)=2.*en*x(l-1)-x(l-2)
      end do
      alamda=0.
      do l=1,nodr
        alamda=alamda+coef(l)*x(l)
      end do
       wlmul(k,j)=wt*(offset+alamda)/(1.+voc)
  171 continue
      go to 7777

c==== legendre case
  702 continue
      read(string(nt1+1:nt2-1),*) iord,iduma,idumb,crval1,cdelt1,
     &  npdum,voc,apb1,apb2,wt,offset,iftype,nodr,pmin,pmax,
     & (coef(l),l=1,nodr)
      do 172 k=1,nump(j)
      p=(k-ltv1)/ltm1_1
      en=(p-(pmax+pmin)/2.)/((pmax-pmin)/2.)
      x(1)=1.
      x(2)=en
      do l=3,nodr
      x(l)=((2*float(l)-3.)*en*x(l-1)-float(l-2)*x(l-2))/float(l-1)
      end do
      alamda=0.
      do l=1,nodr
        alamda=alamda+coef(l)*x(l)
      end do
       wlmul(k,j)=wt*(offset+alamda)/(1.+voc)
  172 continue
      go to 7777

c==== linear spline case
  703 continue
      read(string(nt1+1:nt2-1),*) iord,iduma,idumb,crval1,cdelt1,
     &  npdum,voc,apb1,apb2,wt,offset,iftype,npiece,pmin,pmax,
     & (coef(l),l=1,npiece+1)
      do 173 k=1,nump(j)
      p=(k-ltv1)/ltm1_1
      es=(p-pmin)/(pmax-pmin)*float(npiece)
      jj=int(es)
      aa=(jj+1)-es
      bb=es-jj
      x0=aa
      x1=bb
      alamda=coef(0+jj)*x0+coef(1+jj)*x1
      wlmul(k,j)=wt*(offset+alamda)/(1.+voc)
  173 continue
      go to 7777

c==== cubic spline case
  704 continue
      read(string(nt1+1:nt2-1),*) iord,iduma,idumb,crval1,cdelt1,
     &  npdum,voc,apb1,apb2,wt,offset,iftype,npiece,pmin,pmax,
     & (coef(l),l=1,npiece+3)
      do 174 k=1,nump(j)
      p=(k-ltv1)/ltm1_1
      es=(p-pmin)/(pmax-pmin)*float(npiece)
      jj=int(es)
      aa=(jj+1)-es
      bb=es-jj
      x0=aa**3
      x1=1.+3.*aa*(1.+aa*bb)
      x2=1.+3.*bb*(1.+aa*bb)
      x3=bb**3
      alamda=coef(0+jj)*x0+coef(1+jj)*x1+coef(2+jj)*x2+coef(3+jj)*x3
      wlmul(k,j)=wt*(offset+alamda)/(1.+voc)
  174 continue
      go to 7777

cc==== pixel aray case
c  705 continue
c      go to 7777
cc==== sampled array case
c  706 continue
c      go to 7777

c----------------------------------------------------------------------
 7777 continue
        if(iordp.ne.ipos) then
          if(.NOT.LOREFR) then
            go to 8001
          else
            write(*,*) 'REFR mode: iordp=',iordp,'  ipos=',ipos
            stop
          end if
        end if
 8000 continue
 8001 IF(.NOT.LOREFR) numord=ipos-1

c     check for the 1D case with reduced dimentionality 
      if(nwcsdim.gt.ndim.and.ndim.eq.1) then
        if(iwaxmap01.eq.0) then
          write(*,*) 'why not exist WAXMAP01 despite wcsdim > ndim?!',
     &    nwcsdim,ndim
          stop
        end if 
        read(waxmap01,*) (npln(m),noff(m),m=1,nwcsdim)
        do m=1,nwcsdim
          if(npln(m).eq.0) go to 876
        end do
  876   jsel=noff(m)+1
        npoint=nump(jsel)
        allocate(wla(npoint))
        allocate(spec(npoint))
        do ip=1,npoint
          wla(ip)=wlmul(ip,jsel)
        end do
        npt=npoint
        go to 1003
      end if

c     check the value of norder
      if((.not.lolimit).and.(norder.lt.1.or.norder.gt.numord)) then
c       set fucus on window 0
ccccccccc        result=focusqq(0)
   42 continue
c   42 write(*,'(''Multispe fits file with '',i2,'' orders'')') numord
c      write(*,'(''You can specify the spectrum range to display:'')')
c      write(*,980)(j,int(wlmul(1,j)),int(wlmul(nump(j),j)),j=1,numord)
c  980 format(5(1x,'(',i2.2,':',2i5,')'))
c      write(*,
c     & '(''(Order No.) or (wlmin wlmax) [just enter for all] --> '',$)')
c      read(*,'(a)') linetemp
c      call lablen(linetemp,80,nns,nne)
c      if(nns.gt.nne) then
       wmne=min(wlmul(1,1),wlmul(1,numord))
       wmxe=max(wlmul(nump(1),1),wlmul(nump(numord),numord))
       lolimit=.true.
       go to 652
c      end if
      read(linetemp,*,err=42) norder
      if(norder.lt.100) then
       if(norder.lt.1.or.norder.gt.numord) go to 42
      else
       read(linetemp,*,err=42) wmne,wmxe
       if( wmne.lt.min(wlmul(1,1),wlmul(1,numord))
     &  .or.wmxe.gt.max(wlmul(nump(1),1),wlmul(nump(numord),numord)) )
     &  go to 42
       lolimit=.true.
       norder=0
      end if
  652 continue
      end if
      if((lolimit).and.(norder.lt.1.or.norder.gt.numord)) then
      call bestsel(wmne,wmxe)
      norder=nobest
      end if

c     set norder to nofoc
      IF(.NOT.LOREFR) nofoc=norder
c
      npoint=naxes(1)
c      write(*,*) 'npoint=',npoint
      allocate(wla(npoint))
      allocate(spec(npoint))
      IF(.NOT.LOREFR) THEN
      do k=1,npoint
      wla(k)=wlmul(k,norder)
      end do
      ELSE
      do k=1,npoint
      wla(k)=(w0+(k-1)*dw)/(1.+voc)
      end do
      END IF
      npt=npoint
      go to 1003
 1002 continue
ccc   OTHERS
      write(*,'(''CTYPE1:'',a8)') ctype1
      write(*,*) 'ctype1 unmatch! I cannot read this fits file!'
      stop
      return
 1003 continue

C     initialize variables
c      if(norder.lt.0) npixels=naxes(1)
c      if(norder.eq.0) npixels=naxes(1)
c      if(norder.gt.0) npixels=naxes(1)*naxes(2)
      if(ndim.eq.1) then 
        npixels=naxes(1)
      else 
        npixels=naxes(1)*naxes(2)
      end if
c      allocate(buffer(npixels))

      group=1
      firstpix=1
      nullval=-999
      ncount=0

c      write(*,*) 'npixels=',npixels
      do while (npixels .gt. 0)
C         read all pixels at a time
           nbuffer=npixels
c 6        call ftgpve(unit,group,firstpix,nbuffer,nullval,
c     &            buffer,anynull,status)
  6        call ftgpfe(unit,group,firstpix,nbuffer,
     &            buffer,flagvals,anynull,status)

C         increment pointers and loop back to read the next group of pixels
          npixels=npixels-nbuffer
          firstpix=firstpix+nbuffer
      end do
      if(norder.le.0) then
c     for one-D spec
        do k=1,npt
          spec(k)=buffer(k)
        end do
      else
c     for two-D spec
c
c      first keep all the order data
      IF(.NOT.LOREFR) THEN
      do j=1,numord
        nskip=(j-1)*npt
        do k=1,npt
          spmul(k,j)=buffer(nskip+k)
        end do
      end do
      END IF
c      now set the spectrum for the specified "norder"
        nskip=(norder-1)*npt
        do k=1,npt
          spec(k)=buffer(nskip+k)
        end do
      end if

c_dummyorder
      IF(.NOT.LOREFR.and.norder.le.0) THEN
c     setting nordummy and nump
      if(npt.le.npx) then
        nordummy=1
        nump(1)=npt
      else
        nordummy=npt/(npx-nmar)
        nordummy=nordummy+1
        do j=1,nordummy
          nump(j)=npx
          if(j.eq.nordummy) nump(j)=npt-(npx-nmar)*(nordummy-1)
        end do
      end if
      do j=1,nordummy
      do k=1,nump(j)
      nc=(npx-nmar)*(j-1)+k
      wlmul(k,j)=wla(nc)
      spmul(k,j)=spec(nc)
      end do
      end do
      numord=nordummy
      if(numord.eq.1) norder=1
      if((.not.lolimit).and.(norder.lt.1.or.norder.gt.numord)) then
c  542 write(*,'(''Dummy multispe file with '',i2,'' orders'')') numord
c      write(*,'(''You can specify the spectrum range to display:'')')
c      write(*,5980)(j,int(wlmul(1,j)),int(wlmul(nump(j),j)),j=1,numord)
 5980 format(5(1x,'(',i2.2,':',2i5,')'))
c      write(*,
c     & '(''(Order No.) or (wlmin wlmax) [just enter for all] --> '',$)')
c      read(*,'(a)') linetemp
c      call lablen(linetemp,80,nns,nne)
c      if(nns.gt.nne) then
       wmne=min(wlmul(1,1),wlmul(1,numord))
       wmxe=max(wlmul(nump(1),1),wlmul(nump(numord),numord))
       lolimit=.true.
c       go to 5652
c      end if
c      read(linetemp,*,err=542) norder
c      if(norder.lt.100) then
c       if(norder.lt.1.or.norder.gt.numord) go to 542
c      else
c       read(linetemp,*,err=542) wmne,wmxe
c       if( wmne.lt.min(wlmul(1,1),wlmul(1,numord))
c     &  .or.wmxe.gt.max(wlmul(nump(1),1),wlmul(nump(numord),numord)) )
c     &  go to 542
c       lolimit=.true.
c       norder=0
c      end if
c 5652 continue
      end if
      if((lolimit).and.(norder.lt.1.or.norder.gt.numord)) then
      call bestsel(wmne,wmxe)
      norder=nobest
      end if
c     set norder to nofoc
      nofoc=norder
      END IF
c_dummyorder

ccC     select out spectrum data
ccc
cc      nct=0
cc      do k=1,npt
cc      if((.NOT.LOREFR).AND.(wla(k).lt.wmne)) cycle
cc      if((.NOT.LOREFR).AND.(wla(k).gt.wmxe)) exit
cc      nct=nct+1
cc          if(nct.gt.md) then
cc            if(LOREFR) then
cc              nct=md
cc              go to 79
cc            end if
ccc     too large obspec points  focused order assumed
cc            nptsel=nump(nofoc)
cc            do kk=1,nptsel
cc              wlasel(kk)=wlmul(kk,nofoc)
cc              specsel(kk)=spmul(kk,nofoc)
cc            end do
cc            go to 790
cc          end if
cc      wlasel(nct)=wla(k)
cc      specsel(nct)=spec(k)
cc      end do
cc   79 continue
cc      nptsel=nct
cc      call bestsel(wmne,wmxe)
cc  790 continue

c
C     close the file and free the unit number
 7    call ftclos(unit, status)
      call ftfiou(unit, status)
c     Now deallocating wla,spec,buffer
      deallocate(wla)
      deallocate(spec)
c      deallocate(buffer)
      return
      end

      subroutine lablen(label,nlen,ns,ne)
      implicit real*8 (a-h,o-z)
      character*(*) label
      do n=nlen,1,-1
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) go to 10
      end do
      n=1
  10  ne=n
      do n=1,nlen
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) go to 20
      end do
      n=nlen
  20  ns=n
      return
      end

      subroutine bestsel(wmin,wmax)
      implicit real*8 (a-h,o-z)
      parameter (md=9999)
      common/orderinfo/numord,nump(99),wlmul(md,99),spmul(md,99),
     & nofoc,nobest,novtot,novis(99),xfrac(99)
      data wminsv,wmaxsv/1.d90,1.d90/
      if(wmin.eq.wminsv.and.wmax.eq.wmaxsv) return
      wminsv=wmin
      wmaxsv=wmax
      if(numord.le.1) then
      novtot=numord
      nobest=numord
      nofoc=numord
      return
      end if
      wcen=0.5*(wmin+wmax)
      novtot=0
      nobest=0
      sepmin=1.e30
      do 100 i=1,numord
      if(wmax.lt.wlmul(1,i).or.wmin.gt.wlmul(nump(i),i)) go to 100
      novtot=novtot+1
      novis(novtot)=i
      xfrac(novtot)=(wcen-wlmul(1,i))/(wlmul(nump(i),i)-wlmul(1,i))
      septmp=abs(xfrac(novtot)-0.5)
       if(septmp.lt.sepmin) then
       nobest=novis(novtot)
       sepmin=septmp
      end if
  100 continue
c
      nofoc=nobest
      nofoc=max(1,nofoc)
      nofoc=min(numord,nofoc)
      return
      end
