!___________________________________________________________________________    
!       HI  SURVEYS  SOFTWARE  PACKAGE                                          
!___________________________________________________________________________    
!                                                                               
! This file contains FORTRAN software for use in extracting spectra             
! from the HI surveys of Heiles and Habing (HH), Parkes (PK),                   
! Argentina (AR), Bell Labs (BL), and Weaver and Williams (WW).                 
!                                                                               
! By using these programs and the data that they access, the user               
! agrees to the following:                                                      
! 1. If data from any survey are used in a publication, the user                
!    will reference the original authors of the survey, as if the               
!    original authors had themselves provided the data to the user,             
!      by quoting the appropriate publication in the usual fashion.             
! 2. Whenever data obtained from the NSSDC are used in a publication            
!      acknowledgment shall be made to Dr. Wayne H. Warren Jr. and/or           
!      the Astronomical Data Center/National Space Science Data Center          
!      for storing and distributing the data.                                   
! 3. If the survey data or software described herein are used,                  
!      acknowledgment shall be made to William T. Reach for creating            
!      the tape and writing the software.                                       
!                                                                               
! There is no implied guarantee that these routines will perform the            
! tasks they claim. In other words, user beware--check the results              
! before using them.                                                            
!                                                                               
! The programs included are, in order of appearance:                            
!  DUMP_NSDC, which copies data from tape to disk,                              
!  SPECTRUM_NSDC, which extracts individual spectra from the disk file,         
!  and the subroutine packages                                                  
!  DECODE_NSDC, which contains routines to decode the spectra, and              
!  TAPE HANDLING ROUTINES.                                                      
!                                                                               
! Note: The tape handling routines will not work on any machine other than      
!       a VAX under the VMS operating system; they are included only for        
!       completeness.                                                           
!                                                                               
        program dump_nsdc                                                       
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
!                                                                               
! use this program to copy data from tape to disk                               
! the tape will be positioned to the beginning of the desired survey,           
! based on the assumption that the order of datasets on the tape                
! is HH, PK, AR, BL (for the "high-latitude tape")                              
! (if your tape is ordered differently, change the tape positioning section)    
!                                                                               
        character block*16384,record*4096                                       
        integer spectrum(512)                                                   
        integer nspectra(5),irecord_length(5),nchannels(5)                      
        data nspectra       /  10,  10,  20,  20,  10/                          
        data irecord_length / 643, 420, 396, 703,1204/                          
        data nchannels      / 100,  64,  56, 124, 238/                          
100     write (*,102)                                                           
102     format (/,' NSSDC HI SURVEY TAPE DUMPER',/,                             
     .          x,' (1) Heiles-Habing',/,                                       
     .          x,' (2) Parkes',/,                                              
     .          x,' (3) Argentina',/,                                           
     .          x,' (4) Bell Labs',/,                                           
     .          x,' (5) Weaver-Williams',/,                                     
     .          x,' Enter desired survey: ',$)                                  
        read *,isurvey                                                          
        write (*,'(x,a,$)') 'Which tape drive: '                                
        read *,idrive                                                           
        call getchannel (ichan,idrive)                                          
! position tape at beginning of desired survey                                  
        write (*,'(x,a,$)') 'Does the tape need to be positioned ("y"/n) ? '    
        read (*,'(a)') yn                                                       
        if (yn.eq.'n' .or. yn.eq.'N') goto 200                                  
        if (isurvey.ge.1 .and. isurvey.le.4) then                               
        print *,'Tape labelled "high-latitude" should be mounted'               
        else if (isurvey.eq.5) then                                             
        print *,'Tape labelled "Weaver-Williams" should be mounted'             
        else                                                                    
        print *,'invalid survey...try again'                                    
        goto 100                                                                
        endif                                                                   
        call rewind_tape (ichan)                                                
        if (isurvey.gt.1 .and. isurvey.lt.5) call skipfile(ichan,isurvey-1)     
! position tape at first desired block                                          
! (obtain information from the table of galactic latitude vs. block number)     
200     print *,'first block containing desired data:'                          
        read *,iblock1                                                          
        print *,'number of blocks to read (0=all):'                             
        read *,nblock_out                                                       
        if (iblock1.ne.1) call skipblock (ichan,iblock1-1)                      
        nblock = iblock1-1                                                      
        write (*,250)                                                           
250     format (/,' writing to output file HI; 1 spectrum per record...',/)     
! copy the data from tape to disk                                               
        iok = 0                                                                 
        nspec = 0                                                               
        ib = 0                                                                  
        do while (iok.eq.0)                                                     
                ib = ib + 1                                                     
                if (nblock_out.ne.0 .and. ib.gt.nblock_out) goto 900            
                nbreq = 16384                                                   
                nbytein = 0                                                     
                call readblock (ichan,block,nbreq,nbytein)                      
                if (nbytein.eq.0) then                                          
                        print *,'end of file detected...'                       
                        goto 900                                                
                endif                                                           
                nblock = nblock + 1                                             
                print *,'read block ',nblock,' # bytes = ',nbytein              
                if (ib.eq.1) then                                               
                        irecl = irecord_length(isurvey)                         
                        open (unit=1,file='hi',status='new',recordtype='fixed', 
     .                                  form='formatted',recl=irecl)            
                endif                                                           
                ns = nbytein /irecl                                             
                do i = 1,ns                                                     
                        record(1:irecl) = block((i-1)*irecl+1:i*irecl)          
                        write (1,'(a)') record(1:irecl)                         
                        nspec = nspec + 1                                       
                enddo                                                           
                if (mod(ib,10).eq.0) print *,'finished block ',ib               
        enddo                                                                   
        close (1)                                                               
900     print *,'done...wrote ',nspec,' spectra'                                
        end                                                                     
        program spectrum_nsdc                                                   
c copyright (c) 1989 by William T. Reach, University of California, Berkeley    
!                                                                               
! This routine searches through the HI survey file on disk. It finds the        
! observation nearest a desired position, to within DIST_MAX degrees.           
! The spectrum at this position is written to a three-column file               
! containing the channel number, VLSR of the channel, and antenna               
! temperature, respectively.                                                    
!                                                                               
! The input file is assumed to be assigned to logical name HI; this file        
! should have been produced by the DUMP_NSDC program.                           
!                                                                               
! The output file is the default unit 1 file, FOR001                            
!                                                                               
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer nspectra(5),irecord_length(5),nchannels(5)                      
        data nspectra       /  10,  10,  20,  20,  10/                          
        data irecord_length / 643, 420, 396, 703,1204/                          
        data nchannels      / 100,  64,  56, 124, 238/                          
        real dist_max(5)                                                        
        data dist_max /0.5,1.0,1.0,1.1,0.4/                                     
100     write (*,102)                                                           
102     format (/,' NSSDC HI SURVEY SPECTRUM FINDER',/,                         
     .          x,' (1) Heiles-Habing',/,                                       
     .          x,' (2) Parkes',/,                                              
     .          x,' (3) Argentina',/,                                           
     .          x,' (4) Bell Labs',/,                                           
     .          x,' (5) Weaver-Williams',/,                                     
     .          x,' Enter desired survey: ',$)                                  
        read *,isurvey                                                          
        irecl=irecord_length(isurvey)                                           
        open (unit=20,file='hi',status='old',readonly,                          
     .                  form='formatted',recordtype='fixed',recl=irecl)         
                                                                                
        print *,'desired position (gl,gb):'                                     
        read *,gl0,gb0                                                          
        ifound_it = 0                                                           
        nrec = 0                                                                
        do while (ifound_it.eq.0)                                               
                nrec = nrec + 1                                                 
                read (20,'(a)',end=800) record(1:irecl)                         
                call decode_spectrum (record,isurvey,gl,gb,ra,dec,              
     .                                  vlsr,spectrum)                          
                dist = sqrt((gl-gl0)**2 + (gb-gb0)**2)                          
                if (dist.lt.dist_max(isurvey)) then                             
                        ifound_it = 1                                           
                        print *,'I found it...'                                 
                        print *,'gl = ',gl,'   gb = ',gb                        
                        do i = 1,nchannels(isurvey)                             
                                write (1,210) i,vlsr(i),spectrum(i)             
                        enddo                                                   
                endif                                                           
        enddo                                                                   
        stop                                                                    
210     format (x,i4,2(x,g))                                                    
800     print *,'couldn''t find it...'                                          
        end                                                                     
!___________________________________________________________________________    
!               DECODING  ROUTINES                                              
!                                                                               
! The following routines decode individual spectra of the various HI surveys,   
! returning the antenna temperature as a function of VLSR.                      
!                                                                               
! The first routine, DECODE_SPECTRUM, is a driver for the others and is the     
! only one called by higher-level routines.                                     
!                                                                               
! INPUT VARIABLES:                                                              
! record: character buffer containing the record from the disk file             
! isurvey: the survey number, 1=HH, 2=PK, 3=AR, 4=BL, 5=WW                      
!                                                                               
! OUTPUT VARIABLES:                                                             
! gl,gb: galactic coordinates of the observed position (decimal degrees)        
! ra,dec: equatorial coordinates, decimal hours, degrees                        
! vlsr: array containing the VLSR (km/s) of each spectral channel               
! spectrum: array containing the antenna temperature of each spectral channel   
!                                                                               
        subroutine decode_spectrum (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)  
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer nspectra(5),irecord_length(5),nchannels(5)                      
        data nspectra       /  10,  10,  20,  20,  10/                          
        data irecord_length / 643, 420, 396, 703,1204/                          
        data nchannels      / 100,  64,  56, 124, 238/                          
        if (isurvey.eq.1) then                                                  
                call decode_hh (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)      
        else if (isurvey.eq.2) then                                             
                call decode_pk (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)      
        else if (isurvey.eq.3) then                                             
                call decode_ar (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)      
        else if (isurvey.eq.4) then                                             
                call decode_bl (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)      
        else if (isurvey.eq.5) then                                             
                call decode_ww (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)      
        endif                                                                   
        end                                                                     
!                                                                               
        subroutine decode_hh (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)        
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer*4 ispec(100),igl,igb                                            
        read (record(1:643),901) igl,igb,                                       
     .                          irah,iram,iras,idecd,idecm,idecs,               
     .                          idatey,idatem,idated,ilst,icvel,istat,          
     .                          ispec                                           
901     format (2i6,3i2,i3,2i2,3i2,2i5,i1,1x,100i6)                             
        gl = real(igl)/1000.                                                    
        gb = real(igb)/1000.                                                    
        ra = real(irah) + real(iram)/60. + real(iras)/3600.                     
        dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600.            
        if (idecd.lt.0) dec = -dec                                              
        date = real(idatey) + real(idatem-1)/12. + real(idated-1)/(30.6*12.)    
        cvel = real(icvel)/1000.                                                
        do i = 1,100                                                            
                spectrum(i) = real(ispec(i))/100.                               
        enddo                                                                   
        vwide = (30.e3/1420.4058e6)*2.998e5                                     
        vnarr = (10.e3/1420.4058e6)*2.998e5                                     
        if (cvel.eq.0.) then                                                    
                do i = 1,8                                                      
                        vlsr(i) = real(i-1)*vwide - 90.76                       
                enddo                                                           
                do i = 9,95                                                     
                        vlsr(i) = real(i-9)*vnarr/2. - 45.38                    
                enddo                                                           
                do i = 96,100                                                   
                        vlsr(i) = real(i-96)*vwide + 46.43                      
                enddo                                                           
        else if (cvel.eq.-8.44) then                                            
                do i = 1,8                                                      
                        vlsr(i) = real(i-1)*vwide - 82.32                       
                enddo                                                           
                do i = 9,95                                                     
                        vlsr(i) = real(i-9)*vnarr/2. - 36.94                    
                enddo                                                           
                do i = 96,100                                                   
                        vlsr(i) = real(i-96)*vwide + 54.87                      
                enddo                                                           
        else                                                                    
                print *,'DECODE_HH> unrecognized central velocity'              
        endif                                                                   
        end                                                                     
!                                                                               
        subroutine decode_pk (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)        
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer*4 ispec(64),igl,igb,icv                                         
        read (record(1:420),901) igl,igb,                                       
     .                          irah,iram,iras,idecd,idecm,idecs,               
     .                          icv,istat,                                      
     .                          ispec                                           
901     format (2i6,3i2,i3,2i2,i7,i3,1x,64i6)                                   
        gl = real(igl)/1000.                                                    
        gb = real(igb)/1000.                                                    
        ra = real(irah) + real(iram)/60. + real(iras)/3600.                     
        dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600.            
        if (idecd.lt.0) dec = -dec                                              
        cvel = real(icv)/1000.                                                  
        do i = 1,64                                                             
        spectrum(i) = real(ispec(i))/100.                                       
        vlsr(i) = cvel + real(i-1)*6.99                                         
        enddo                                                                   
        end                                                                     
!                                                                               
        subroutine decode_ar (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)        
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer*4 ispec(56),ibadscan(5),igl,igb,iv1,iv56                        
        read (record(1:396),901) igl,igb,                                       
     .                          irah,iram,iras,idecd,idecm,idecs,               
     .                          ilsth,ilstm,ilsts,idatey,idatem,idated,         
     .                          (ibadscan(k),k=1,5),                            
     .                          iv1,iv56,istat,                                 
     .                          ispec                                           
901     format (2i6,3i2,i3,2i2,6i2,5i2,2i6,i1,56i6)                             
        gl = real(igl)/1000.                                                    
        gb = real(igb)/1000.                                                    
        ra = real(irah) + real(iram)/60. + real(iras)/3600.                     
        dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600.            
        if (idecd.lt.0) dec = -dec                                              
        date = real(idatey) + real(idatem-1)/12. + real(idated-1)/(30.6*12.)    
        v1 = real(iv1)/100.                                                     
        v56 = real(iv56)/100.                                                   
        dv = (v56-v1)/55.                                                       
        do i = 1,56                                                             
                spectrum(i) = real(ispec(i))/100.                               
                vlsr(i) = v1 + (i-1)*dv                                         
        enddo                                                                   
        end                                                                     
!                                                                               
        subroutine decode_bl (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)        
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer*4 ispec(124),igl,igb,icv                                        
        integer*4 icoeffn(3),icoeffw(3)                                         
        real coeffn(3),coeffw(3)                                                
        read (record(1:703),901) igl,igb,                                       
     .                          irah,iram,iras,idecd,idecm,idecs,               
     .                          ich0,icv,itint,irms,nchbase,                    
     .                          icoeffn,icoeffw,                                
     .                          ispec                                           
901     format (2i6,3i2,i3,2i2,i4,i6,i4,i5,i3,6i6,124i5)                        
        gl = real(igl)/100.                                                     
        gb = real(igb)/100.                                                     
        ra = real(irah) + real(iram)/60. + real(iras)/3600.                     
        dec = abs(real(idecd)) + real(idecm)/60. + real(idecs)/3600.            
        if (idecd.lt.0) dec = -dec                                              
        ch0 = real(ich0)/10.                                                    
        cvel = real(icv)/100.                                                   
        do i = 1,3                                                              
                coeffn(i) = real(icoeffn(i))/1.e4                               
                coeffw(i) = real(icoeffw(i))/1.e4                               
                enddo                                                           
                do i = 1,124                                                    
                        spectrum(i) = real(ispec(i))/100.                       
                        vlsr(i) = cvel + (real(i)-ch0)*5.2767                   
        enddo                                                                   
        end                                                                     
!                                                                               
        subroutine decode_ww (record,isurvey,gl,gb,ra,dec,vlsr,spectrum)        
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        character record*4096                                                   
        real spectrum(512),vlsr(512)                                            
        integer*4 ispec(238),icv,igl,igb                                        
        read (record(1:1204),901) igl,igb,icv,ispec                             
901     format (i4,i5,i4,i1,238i5)                                              
        gl = real(igl)/10.                                                      
        gb = real(igb)/100.                                                     
        cvel = real(icv)/10.                                                    
        do i = 1,238                                                            
                spectrum(i) = real(ispec(i))/10.                                
                vlsr(i) = cvel + (real(i)-122)*1.0553                           
        enddo                                                                   
        return                                                                  
        do i = 1,238                                                            
                spectrum(i) = 0.                                                
                vlsr(i) = 0.                                                    
        enddo                                                                   
        end                                                                     
!______________________________________________________________________________ 
!               TAPE  HANDLING  ROUTINES                                        
!                                                                               
! The following routines are a tape driver package that perform the interface   
! between FORTRAN and the VAX/VMS system I/O routines. They will not work on    
! any other operating system. They assume the tape drives are assigned to       
! the names MUB0: and MUB1:, and that there are two tape drives. If your        
! system is configured differently, change the routine GETCHANNEL accordingly.  
!                                                                               
c*****************************************************************GETCHANNEL    
c                                                                               
c ASSIGN CHANNEL TO THE TAPE DRIVE                                              
c N.B. ALL I/O OPERATIONS TO TAPE DRIVE REQUIRE A CHANNEL                       
c                                                                               
        subroutine getchannel (channel,drive)                                   
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        integer*4 status,sys$assign,channel,drive                               
        if (drive .eq. -1) then                                                 
          status = sys$assign ('MUB0:',channel,,)                               
          if (.not. status) status = sys$assign ('MUB1:',channel,,)             
          if (.not. status) call lib$signal(%val(status))                       
          return                                                                
         endif                                                                  
        if (drive .gt. 2 .or. drive .lt. -1) then                               
          print *,'tape drive ',drive,' does not exist .. access denied'        
          return                                                                
         endif                                                                  
        if (drive .eq. 0) status = sys$assign ('MUB0:',channel,,)               
        if (drive .eq. 1) status = sys$assign ('MUB1:',channel,,)               
        if (.not. status) call lib$signal(%val(status))                         
        end                                                                     
c                                                                               
c******************************************************************READBLOCK    
c                                                                               
c READ NEXT LOGICAL BLOCK ON TAPE                                               
c BUFFER = CHARACTER VARIABLE TO RECIEVE BLOCK                                  
c SIZE = NUMBER OF BYTES TO RECEIVE                                             
c NBYTES = NUMBER OF BYTES ACTUALLY READ INTO BUFFER                            
c                                                                               
        subroutine readblock (channel,buffer,size,nbytes)                       
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        integer*4 status,sys$qiow,iosb(2),channel,size                          
        integer*2 func                                                          
        character*(*) buffer                                                    
        func = 33                                                               
        status = sys$qiow(, %val(channel), %val(func), iosb, , ,                
     +   %ref(buffer(1:1)), %val(size), , , , )                                 
        if (.not. status) call lib$signal(%val(status))                         
        do 10 n = 30,12,-1                                                      
          if (iosb(1) .lt. 2**n) goto 10                                        
          iosb(1) = iosb(1) - 2**n                                              
          nbytes = nbytes + 2**(n-16)                                           
10      continue                                                                
        return                                                                  
        end                                                                     
c                                                                               
c*****************************************************************REWIND_TAPE   
c                                                                               
c REWIND TAPE                                                                   
c                                                                               
        subroutine rewind_tape (channel)                                        
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        integer*4 status,sys$qiow,iosb(2),channel                               
        integer*2 func                                                          
        func = 36                                                               
        status = sys$qiow(,%val(channel),%val(func),iosb,,,,,,,,)               
        if (.not. status) call lib$signal(%val(status))                         
        return                                                                  
        end                                                                     
c                                                                               
c*****************************************************************SKIPBLOCK     
c                                                                               
        subroutine skipblock (channel,nskip)                                    
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        integer*4 status,sys$assign,sys$qiow,iosb(2),channel                    
        integer func*2,nskip*2                                                  
        func = 38                                                               
        status = sys$qiow (,%val(channel),%val(func),iosb,,,                    
     +                          %val(nskip),%val(1),,,,)                        
        if (.not. status) call lib$signal(%val(status))                         
        return                                                                  
        end                                                                     
c                                                                               
c*****************************************************************SKIPFILE      
c                                                                               
        subroutine skipfile (channel,nskip)                                     
! copyright (c) 1989 by William T. Reach, University of California, Berkeley    
        integer*4 status,sys$assign,sys$qiow,iosb(2),channel,nskip              
        integer func*2                                                          
        func = 37                                                               
        status = sys$qiow (,%val(channel),%val(func),iosb,,,                    
     +                          %val(nskip),%val(1),,,,)                        
        if (.not. status) call lib$signal(%val(status))                         
        return                                                                  
        end                                                                     
