Implicit None
	Character*32 file1,file2*10,file3,file4,header(200),outfile*11
	character *20 junk
	Integer i,j,z,num,nsteps,nnod,year,month,day,nx,nz,ny
	integer posh,posv,aa,bb,cc,ii,jj,y,l,zzz,maxlatdeg,bbb
	integer aainc,bbinc,count(6000),latbins,longbins
	integer lz,il,pixels
	Real  scane,scann,offsete,offsetn,b,m
	real  spec(3000,6000),delta,wave(6000)
	real  currentspec(6000)
	real  waveset(6000)
	real*8  yfact,bfact,fact,R,newwave(6000),minw,maxw
	real   maxwave,minwave,pix,xxx
	
***    spec(nsteps*ny,nx)
***    wave(nx)
**     newwave 
***     newspec(1) is the counter to tell how many spec have been added
***     aainc and bbinc are the numbers used a increments in final map long and lat
***     I have modified this so that I throw away the top and bottom two rows of
***	data along the slit( they are noisier than the stuff inthe middle)
****	Read in data file holding script

	R=600000.		!resolving power
	print*, 'Please input the name of the script file!'
	read*,file1
	open(10,file=file1,status='old')
	read(10,*) num

****** opens this file to create the wavelength grid
****** need to put the first file in script twice (once to give wavegrid 2nd for data use)
****** Only put in num of different files in script since num is used after setting up
****** the wavenumber grid.
		read(10,'(1a9)') file2
        	file3=file2(1:9)//'.hd'
        	open(30,file=file3,status='old')

        	do i=1,200
       		read(30,'(1a32)',end=80) header(i)
        	end do
   80   read(header(8)(11:14),*) scane
        read(header(8)(20:23),*) scann
        read(header(9)(11:14),*) offsete
        read(header(9)(20:23),*) offsetn
        read(header(11)(11:12),*) nsteps
        read(header(18)(11:12),*) nnod
        read(header(45)(11:18),'(i2,1x,i2,1x,i2)') year,month,day
        read(header(i-4)(13:16),*) nx
        read(header(i-3)(15:16),*) ny
        read(header(i-2)(11:11),*) nz
        close(30)

        open(60,file=file2,access='direct',recl=nx*4,
     &                  form='unformatted')
        read(60,rec=1) (waveset(y),y=1,nx)
        close(60)
	

        maxw=waveset(nx)+.7
        minw=waveset(1)-.7
        delta=minw/R
	pixels=1+((maxw-minw)/delta)
		print*,pixels
        newwave(1)=minw
        do ii=2,pixels
           newwave(ii)=(dexp(-1*(dlog(R)-dlog(newwave(ii-1)))))
     &+newwave(ii-1)
        end do


********  start working

	do j=1,num
	print*,'file',j

**** Read in header of file j+1 in list

                read(10,'(1a9)') file2
                file3=file2(1:9)//'.hd'
                open(30,file=file3,status='old')

                do i=1,200
                read(30,'(1a32)',end=90) header(i)
                end do
   90   read(header(8)(11:14),*) scane
        read(header(8)(20:23),*) scann
        read(header(9)(11:14),*) offsete
        read(header(9)(20:23),*) offsetn
        read(header(11)(11:12),*) nsteps
        read(header(18)(11:12),*) nnod
        read(header(45)(11:18),'(i2,1x,i2,1x,i2)') year,month,day
        read(header(i-4)(13:16),*) nx
        read(header(i-3)(15:16),*) ny
        read(header(i-2)(11:11),*) nz
        close(30)

***** Read in all the data file for one scan

	open(60,file=file2,access='direct',recl=nx*4,
     &                  form='unformatted')
        read(60,rec=1) (wave(y),y=1,nx)

	do z=1,ny*nz
        read(60,rec=z) (spec(z,y),y=1,nx)
	end do
        close(60)


	
	junk=file2(1:9)//'.dat'	
        open(70,file=junk,access='direct',recl=(pixels)*4,
     &                  form='unformatted')

	do z=1,ny*nz
	xxx=ny*nz
	if ((((real(z)-1.)/real(ny))-((z-1)/ny)).ne.0) then
*** This laces together each spectrum
                do ii=1,pixels
			count(ii)=0
			currentspec(ii)=0
                        do jj=1,nx-1
                        if (spec(z,jj).ne.0) then
			if (newwave(ii).eq.wave(jj)) then
                        if (count(ii).eq.1) then
                        pix=(jj/256.-jj/256)*256.
                        do il=1,80
				if (((jj-il)/256.-(jj-il)/256).eq.0) then
                                minwave=wave(jj-il+1)
                                goto 395
				end if
                                if (spec(z,jj-il).eq.0) then
                                minwave=wave(jj-il+1)
                                goto 395
                                end if
                        end do
 395                    do lz=1,80
                                if (spec(z,256*(jj/256)-lz).ne.0) then
                                maxwave=wave(256*(jj/256)-lz)
                                goto 396
                                end if
                        end do
 396                    fact=1./(maxwave-minwave)
                        bfact=1.-fact*maxwave
                        yfact=fact*wave(jj)+bfact
                m=(spec(z,jj+1)-spec(z,jj))/
     &(wave(jj+1)-wave(jj))
                        b=spec(z,jj)-m*wave(jj)
                        currentspec(ii)=currentspec(ii)*(1.-yfact)
     &                          +(m*newwave(ii)+b)*(yfact)
                        else
			currentspec(ii)=currentspec(ii)+spec(z,jj)
			count(ii)=count(ii)+1
			end if
			end if
                	if (newwave(ii).gt.wave(jj).and.
     &newwave(ii).lt.wave(jj+1)) then
			if (count(ii).eq.1) then
			pix=(jj/256.-jj/256)*256.
			do il=1,80
				if (((jj-il)/256.-(jj-il)/256).eq.0) then
                                minwave=wave(jj-il+1)
                                goto 393
				end if
				if (spec(z,jj-il).eq.0) then
				minwave=wave(jj-il+1)
				goto 393
				end if
			end do
 393			do lz=1,80
				if (spec(z,256*(jj/256)-lz).ne.0) then
				maxwave=wave(256*(jj/256)-lz)
				goto 394
				end if
			end do
 394			fact=1./(maxwave-minwave)
			bfact=1.-fact*maxwave
			yfact=fact*wave(jj)+bfact	
                m=(spec(z,jj+1)-spec(z,jj))/
     &(wave(jj+1)-wave(jj))
                        b=spec(z,jj)-m*wave(jj)
                        currentspec(ii)=currentspec(ii)*(1.-yfact)
     &                          +(m*newwave(ii)+b)*(yfact)
			else
                m=(spec(z,jj+1)-spec(z,jj))/
     &(wave(jj+1)-wave(jj))
                        b=spec(z,jj)-m*wave(jj)
                        currentspec(ii)=currentspec(ii)
     &            		+(m*newwave(ii)+b)
			count(ii)=count(ii)+1
			end if
                        end if
                        end if
                        end do
                        end do
        write(70,rec=z) (currentspec(zzz),zzz=1,pixels)
			else
	write(70,rec=z) (real(newwave(zzz)),zzz=1,pixels)
			end if
		end do
		close(70)
	end do
  	end