* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      SUBROUTINE NIC_SIMUL(LINE,ERROR)
C ----------------------------------------------------------------------
C 	Support routine for command
C	SIMUL  /scan scan_number /input xxxx
C or    SIMUL  /scan scan_number /GAUSS x y zmax b
C       x, y ,b in arcseconds
C
C Arguments :
C	LINE	C*(*)	Command line			Input
C	ERROR	L	Logical error flag		Output
C ----------------------------------------------------------------------
      CHARACTER*(*) LINE
      LOGICAL ERROR
*
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
*
      INTEGER*4 IPMAP,IADDR
      INTEGER*4 POINTER
      REAL*8          BEAMSEP
*
      LOGICAL   SIC_PRESENT, DATA_TYPE
      CHARACTER CHAIN*80,FICH*80, FTYPE*32, EXTEN*4, STRING*40
      CHARACTER TYPE*4
      PARAMETER   (EXTEN='.cnv')
      INTEGER   NFILE,LENC,FLEN,NC,I,J
      REAL      REF1,REF2,VAL1,VAL2,INC1,INC2,DX,DY
      INTEGER   I_SCAN,I_POINT,NUMBER
      INTEGER   BLC(4),TRC(4)
      INTEGER   SIC_GETVM,IER
      REAL      X,Y,R2,A,Z
      REAL      X0,Y0,ZMAX,B,B2,CSTE
      INTEGER   NX,NY
      LOGICAL   EXIST_NMB,UPDATE_NMB
      INTEGER   SIMUL_TYPE
      REAL*8      ZIGZAG
*
      DATA      BLC /0,0,0,0/, TRC /0,0,0,0/
*
* Check if all parameters are present
*
C      CALL MESSAGE (2,1,'SIMULATE',
C     $'Data simulation  is experimental ...')
      ERROR = .FALSE.
      IF (.NOT.SIC_PRESENT(1,0) .OR. .NOT.SIC_PRESENT(1,1))
     $ERROR = .TRUE.
      IF (.NOT.SIC_PRESENT(2,0) .AND..NOT.SIC_PRESENT(3,0))
     $ERROR = .TRUE.
      IF (SIC_PRESENT(2,0) .AND. .NOT.SIC_PRESENT(2,1))  ERROR = .TRUE.
      IF (SIC_PRESENT(2,0).AND.SIC_PRESENT(3,0) ) THEN
         CALL MESSAGE(6,3,'SIMULATE',
     $   'INPUT and GAUSS options are mutually exclusive')
         ERROR=.TRUE.
         RETURN
      ENDIF
      IF (SIC_PRESENT(3,0) .AND.
     $(.NOT.SIC_PRESENT(3,1) .OR. .NOT.SIC_PRESENT(3,2) .OR.
     $.NOT.SIC_PRESENT(3,3) .OR. .NOT.SIC_PRESENT(3,4))) THEN
         CALL MESSAGE(6,3,'SIMULATE',
     $   'SIMUL  /SCAN scan_number /GAUSS x y zmax b')
         RETURN
      ENDIF
      IF (ERROR) THEN
         CALL MESSAGE(6,3,'SIMULATE','Missing options')
         RETURN
      ENDIF
*
* Read nmb file
*
      CALL SIC_CH (LINE,1,1,CHAIN,NC,.TRUE.,ERROR)
      IF (ERROR) GOTO 99
      IF (INDEX(CHAIN,'.').EQ.0) WRITE (6,'('' Defaulted to .NMB'')')
      CALL SIC_PARSEF(CHAIN,FICH,'BOLO_DATA:',DEFEXT)
*      CALL SIC_PARSEF(CHAIN,FICH,' ',DEFEXT)
      I_SCAN = INDEX(FICH,CHAIN(1:LENC(CHAIN)))
      I_POINT = INDEX(FICH(I_SCAN:),'.')
      FTYPE = FICH(I_SCAN+I_POINT-1:)
      FLEN = LENC(FICH(I_SCAN+I_POINT-1:))
      CALL CHECK_EXTENSION ( FTYPE, FLEN, ERROR)
      IF (ERROR) RETURN
      IF (FTYPE.EQ.'.NMB'.OR.FTYPE.EQ.'.nmb') THEN
         READ (CHAIN,'(i4)',ERR=99) NUMBER
         WRITE (FICH(I_SCAN:),'(i4.4)') NUMBER ! format i4.4
         WRITE(FICH(I_SCAN+4:),'(a4)')'.nmb'
         NFILE = LENC(FICH)
         DATA_TYPE = .FALSE.
         CALL  INPUT_NMB(FICH,NFILE,ERROR,DATA_TYPE)
         INQUIRE (FILE=FICH(1:NFILE),EXIST=EXIST_NMB)
         IF (.NOT.EXIST_NMB) RETURN
      ELSEIF (FTYPE.EQ.'.ASC'.OR.FTYPE.EQ.'.asc') THEN
         CALL  INPUT_ASC(FICH,NFILE,ERROR,DATA_TYPE)
      ENDIF
      IF (ERROR) RETURN
* Option to update signal in nmb file
      UPDATE_NMB = .FALSE.
      IF (SIC_PRESENT(1,2)) THEN
         CALL SIC_KE (LINE,1,2,CHAIN,NC,.TRUE.,ERROR)
         IF (CHAIN(1:2).EQ.'UP') UPDATE_NMB = .TRUE.
      ENDIF
*
* Read map
*
      IF (SIC_PRESENT(2,0)) THEN
         SIMUL_TYPE = 1
         CALL SIC_CH (LINE,2,1,CHAIN,NC,.TRUE.,ERROR)
         IF (ERROR) GOTO 99
         CALL SIC_LOWER(CHAIN)
         I = INDEX(CHAIN,'.')
         IF (I.EQ.0) THEN
            WRITE (6,'('' Defaulted to .cnv'')')
            TYPE = 'cnv'
         ELSE
            TYPE = CHAIN(I+1:)
         ENDIF
         CALL SIC_PARSEF(CHAIN,X_FILE,' ',EXTEN)
         CALL GDF_GEIS  (X_ISLO,ERROR)   ! reserve slot
         CALL GDF_REIS  (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
         IF (ERROR) THEN
            WRITE (STRING,'(a,1x,a)') 'Error on file ',X_FILE(1:NC)
            CALL MESSAGE (2,1,'SIMUL',STRING)
            CALL GDF_FRIS (X_ISLO,ERROR) ! free slot
            RETURN
         ENDIF
         CALL GDF_READX  (X_ISLO,ERROR)
         CALL GDF_GEMS  (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
         IPMAP = POINTER (X_ADDR,MEMORY)
         REF1 = X_REF1
         REF2 = X_REF2
         VAL1 = X_VAL1/SEC_TO_RAD
         VAL2 = X_VAL2/SEC_TO_RAD
         INC1 = X_INC1/SEC_TO_RAD
         INC2 = X_INC2/SEC_TO_RAD
         NX = X_DIM(1)
         NY = X_DIM(2)
      ENDIF
*
* Build gaussian map
*
      IF (SIC_PRESENT(3,0)) THEN
         SIMUL_TYPE = 2
         CALL SIC_R4 (LINE,3,1,X0,.TRUE.,ERROR)
         CALL SIC_R4 (LINE,3,2,Y0,.TRUE.,ERROR)
         CALL SIC_R4 (LINE,3,3,ZMAX,.TRUE.,ERROR)
         CALL SIC_R4 (LINE,3,4,B,.TRUE.,ERROR)
         CALL SIC_GET_REAL('MAP_CELL[1]',INC1,ERROR)
         CALL SIC_GET_REAL('MAP_CELL[2]',INC2,ERROR)
         DX = SQRT(DAZM**2+4*OAZM**2)* 3600    ! in arcsec
         DY = DX
         NX = NINT(DX/INC1)
         NX = 2*(NX/2)+1
         NY = NINT(DY/INC2)
         NY = 2*(NY/2)+1
         VAL1 = 0
         VAL2 = 0
         REF1 = NX/2+1
         REF2 = NY/2+1
         IER = SIC_GETVM(NX*NY,IADDR)    ! gaussian map
         IPMAP = POINTER(IADDR,MEMORY)
         IF (IER.NE.1) THEN
            WRITE(6,*)' *** Unable to allocate memory ***'
         ENDIF
!  model = zmax* exp (-4*ln(2)*(r/b)**2)
         CSTE = -4*LOG(2.)
         B2   = B**2
         OPEN(UNIT=7,FILE ='t.t')
         DO J=1,NY
            WRITE(7,*) J
            Y = VAL2 + (J-REF2)*INC2
            DO I=1,NX
               X = VAL1 + (I-REF1)*INC1
               R2 = (X-X0)**2+(Y-Y0)**2
               A = CSTE*R2/B2
               Z = 0
               IF (A.GT.-20.) Z = ZMAX*EXP(A)
               CALL R4TOR4(Z,MEMORY(IPMAP+I-1+(J-1)*NX),1)
               WRITE(7,'(e8.2,$)') Z
            ENDDO
            WRITE(7,*)' '
         ENDDO
         CLOSE(UNIT=7)
      ENDIF
*  Compute signal
      CALL  COMPUTE_SIGNAL_FROM_MAP(
     $MEMORY(IPMAP),VAL1,INC1,REF1,VAL2,INC2,REF2,
     $NX,NY,
     $MEMORY(IP_LST),MEMORY(IP_EL),MEMORY(IP_SCAN_COORD),
     $MEMORY(IP_SCAN_COORD+NRECORD),MEMORY(IP_SIGNAL),
     $MEMORY(IP_AZOFF),MEMORY(IP_ELOFF),MEMORY(IP_FLAG),
     $MEMORY(IPSUBSCAN+NSUBSCAN),ERROR)
*
*
      IF (UPDATE_NMB.AND..NOT.ERROR) THEN
         IF (SIMUL_TYPE.EQ.2..AND.ZMAX.LT.10) THEN
            WRITE(6,*)
     $      ' Low amplitude : cannot be written',
     $      ' in nmb format'
         ELSE
            WRITE(6,*)' Updating nmb file ',FICH
            WRITE(6,*)
     $      ' Amplitudes must be in the range of bolometer channels ',
     $      ' (from 100 to 200000) to be stored in nmb format'
            CALL SAVE_NMB(FICH,MEMORY(IPSUBSCAN+NSUBSCAN),
     $      MEMORY(IPSUBSCAN+2*NSUBSCAN),
     $      MEMORY(IP_LST),MEMORY(IP_AZ),MEMORY(IP_EL),
     $      MEMORY(IP_SCAN_COORD),
     $      MEMORY(IP_SCAN_COORD+NRECORD),
     $      MEMORY(IP_SIGNAL))
         ENDIF
      ENDIF
      ERROR=.FALSE.
*
*   Free memory
*
      IF (SIMUL_TYPE.EQ.1) THEN
         IF (SIC_PRESENT(2,0)) CALL GDF_FRIS (X_ISLO,ERROR)
      ELSE
         CALL FREE_VM(NX*NY,IADDR)
      ENDIF
      BEAMSEP = WOBBLER_THROW
      CALL SIC_LET_DBLE('BEAMSEP',BEAMSEP,ERROR)
      ZIGZAG = 0
      CALL SIC_LET_DBLE('zigzag',ZIGZAG,ERROR)
      RETURN
*
99    ERROR=.TRUE.
      RETURN
      END
*
C ----------------------------------------------------------------------
C
C     COMPUTE_SIGNAL_FROM_MAP
C
C ----------------------------------------------------------------------
*
      SUBROUTINE  COMPUTE_SIGNAL_FROM_MAP(RADEC_MAP,VAL1,INC1,REF1,
     $VAL2,INC2,REF2,NX,NY,LST,EL,XCOO,YCOO,SIGNAL,
     $OA,OE,FLAG,EPOCH_SUBSCAN,ERROR)
*
      INCLUDE 'const.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
*
      REAL      RADEC_MAP(1)
      REAL      EL(1),LST(1),XCOO(1),YCOO(1),SIGNAL(1)
      REAL      OA(1),OE(1)
      REAL      VAL1,INC1,REF1
      REAL      VAL2,INC2,REF2
      INTEGER   NX,NY
      LOGICAL   FLAG(1),ERROR
      REAL      EPOCH_SUBSCAN(1)
* local variables
      REAL*8    DRA,DDEC,DLST,DAZ,DEL,LATITUDE_DEG
      REAL*8    EPOCHFR,EPOCHTO,AMAT(9),AVEC(3),BVEC(3)
      REAL*8    RA_FROM,DEC_FROM,RA_DAY,DEC_DAY
      INTEGER   I,J,K,ICHAN
      REAL      FUNWID,TOTWIDX,TOTWIDY,POSX,POSY,Z
      INTEGER   IFIRST_X,ILAST_X,IFIRST_Y,ILAST_Y
      REAL      COS_EL,SIN_EL
      REAL      B2,R2,CSTE,A
      REAL      FLUX(2),WEIGHT,VAL
      INTEGER   IRECORD,IWOB,LAST_CHAN
      LOGICAL FOUND, SIC_CTRLC
* Initialization
      ERROR = .FALSE.
      FUNWID = VLIGHT/FREQUENCY /(TELESCOPE_DIAMETER)
     $/SEC_TO_RAD/1.D9                   ! in arcsecs
!                                                ! no magic number ...
      TOTWIDX = ABS(4*FUNWID/INC1)
      TOTWIDY = ABS(4*FUNWID/INC2)
      B2 = ABS(FUNWID*FUNWID/INC1/INC2)
      CSTE = -4*LOG(2.)
* Compute RA_DAY and DEC_DAY
      RA_FROM = SOURCE_LAMBDA
      DEC_FROM = SOURCE_BETA
      LATITUDE_DEG = SITE_LATITUDE/DEG_TO_RAD
      EPOCHFR = EPOCH
      CALL R4TOR8(EPOCH_SUBSCAN(1),EPOCHTO,1)
      CALL PRECMAT(EPOCHFR,EPOCHTO,AMAT) ! compute precession matrix
      CALL DCOSIN(RA_FROM,DEC_FROM,AVEC) ! compute vector
      CALL MATMUL(AMAT,AVEC,BVEC,-9)     ! multiply by precession matrix
      CALL DANGLE(RA_DAY,DEC_DAY,BVEC)   !
* compute signal
c     OPEN(UNIT=7,FILE='w.w')
      DO ICHAN=1,NCHAN
         FOUND = SIC_CTRLC ()
         LAST_CHAN=ICHAN+4
         IF (LAST_CHAN.GT.NCHAN) LAST_CHAN=NCHAN
         IF (MOD(ICHAN,5).EQ.1) WRITE(6,*)
     $   'Computing signal for channels ',ICHAN,
     $   ' to ',LAST_CHAN,' ...'
         IF (.NOT.FLAG(ICHAN)) THEN
            DO IRECORD= 1,NRECORD
               DLST = LST(IRECORD)
               COS_EL = COS(FIELD_ROTATION*EL(IRECORD)*DEG_TO_RAD)
               SIN_EL = SIN(FIELD_ROTATION*EL(IRECORD)*DEG_TO_RAD)
               DO IWOB =1,2
                  DAZ = XCOO(IRECORD)/10.
     $            + (-OA(ICHAN))*COS_EL + (-OE(ICHAN))*SIN_EL
                  DEL  = YCOO(IRECORD)/10.
     $            - (-OA(ICHAN))*SIN_EL + (-OE(ICHAN))*COS_EL
                  IF (IWOB.EQ.1) DAZ = DAZ - WOBBLER_THROW/2.
                  IF (IWOB.EQ.2) DAZ = DAZ + WOBBLER_THROW/2.
                  CALL NIC_AZEL_TO_RADEC(DAZ,DEL,DLST,LATITUDE_DEG,RA_DA
     $            Y,
     $            DEC_DAY,DRA,DDEC)
c            write(6,*)daz,del,dlst,latitude_deg,ra_day,dec_day,
c     $           dra,ddec
                  POSX = (DRA-VAL1) /INC1 + REF1
                  POSY = (DDEC-VAL2)/INC2 + REF2
!            write(6,*) irecord,iwob,daz,del,dra,ddec
                  IFIRST_Y = POSY - TOTWIDY/2.+1.
                  ILAST_Y  = POSY + TOTWIDY/2.
                  IF (IFIRST_Y.LT.1)   IFIRST_Y = 1
                  IF (ILAST_Y.GT.NY) ILAST_Y  = NY
                  IFIRST_X = POSX - TOTWIDX + 1.
                  ILAST_X  = POSX + TOTWIDX
                  IF (IFIRST_X.LT.1)    IFIRST_X = 1
                  IF (ILAST_X .GT.NX) ILAST_X  = NX
                  FLUX(IWOB) = 0
                  WEIGHT = 0
                  DO J = IFIRST_Y, ILAST_Y
                     DO I=IFIRST_X, ILAST_X
                        R2 = (I-POSX)**2+(J-POSY)**2
                        A = CSTE *R2/B2
                        Z = 0
                        IF (A.GT.-20) Z = EXP(A)
                        IF (Z.NE.0) THEN
                           K = I + (J-1)*NX
                           CALL R4TOR4(RADEC_MAP(K-1),VAL,1)
                           IF (ABS(VAL-BLANKING).GT.1) THEN
                              FLUX(IWOB) = FLUX(IWOB) + VAL*Z
                              WEIGHT = WEIGHT+Z
                           ENDIF
                        ENDIF
                     ENDDO
                  ENDDO
                  IF (WEIGHT.GT.1.E-6) THEN
                     FLUX(IWOB) = FLUX(IWOB)/WEIGHT
                  ELSE
                     FLUX(IWOB) = BLANKING
                  ENDIF
c                  IF (ICHAN.EQ.1) THEN
c                     WRITE(7,*) IRECORD,IWOB,FLUX(IWOB),WEIGHT,
c     $               IFIRST_Y, ILAST_Y,IFIRST_X, ILAST_X,POSX,POSY
c                  ENDIF
               ENDDO                     ! loop iwob=1,2
               IF (ABS(FLUX(1)-BLANKING).GT.1.AND.
     $         ABS(FLUX(2)-BLANKING).GT.1) THEN
                  FLUX(1) = FLUX(2)-FLUX(1)
               ELSE
                  FLUX(1) = BLANKING
               ENDIF
               CALL R4TOR4(FLUX(1),SIGNAL(IRECORD+(ICHAN-1)*NRECORD),1)
               IF (FOUND) GOTO 300
            ENDDO                        ! loop irecord
         ENDIF                           ! if .not flag
      ENDDO                              ! loop ichan
c      CLOSE(UNIT=7)
      RETURN
300   ERROR=.TRUE.
      RETURN
      END
*
* ----------------------------------------------------------------------
*
*     NIC_FUNGEN
*
*       this subroutine sets up the convolving function used to         *
*       filter and regrid temporal data.                                *
* ----------------------------------------------------------------------
*
      SUBROUTINE NIC_FUNGEN(FUNC,IFUN)
*
      INCLUDE 'const.inc'
      REAL*8    FUNC(1), TAP, XF0, XF1
      INTEGER   I,IFUN
*
      IF (IFUN.GT.2.OR.IFUN.LT.0) THEN
         WRITE (6,'(a)') 'No tapering function applied'
      ENDIF
      TAP = 1D0
      XF0 = PI/5D1
      XF1 = 0D0
      FUNC(1) = 1D0
      DO I = 2, 201
         XF1 = XF1 + XF0
         TAP = 1D0
*
* ---- Tapering with a linear slope in the last quarter
*
*         if (i.ge.151) then
*            tap = (201d0 - i)/5d1
*         endif
*
* ---- Tapering with a parabolic profile
*
         IF (IFUN.EQ.0) TAP = 1-((I-1)/2D2)**2
*
* ---- Tapering with a gaussian profile
*
*
         IF (IFUN.EQ.1) TAP = EXP(-((I-1)/1D2)**2)
*
* ---- Tapering with a linear slope
*
         IF (IFUN.EQ.2) TAP = 1-((I-1)/2D2)
*
         FUNC(I) = TAP * SIN(XF1)/XF1
      ENDDO
      RETURN
      END
*
*
C******************************************************************
C     NIC_RADEC_TO_AZEL
C
C
C   THIS ROUTINE TAKES AS INPUT THE RA AND DEC OFFSET RELATIVE     *
C   TO THE FIELD CENTRE, THE L.S.T., THE OBSERVATORY LATITUDE     *
C   AND THE CURRENT EPOCH EQUATORIAL COORDINATES OF THE FIELD     *
C   CENTRE AND RETURNS THE AZ AND EL OFFSETS RELATIVE TO THE     *
C   FIELD CENTRE.                                                 *
C   DRA,DDEC  = OFFSETS RELATIVE TO FIELD CENTRE IN RA-DEC
C   LST = L.S.T. OF DATA POINT                                   *
C   SITE_LATITUDE = OBSERVATORY LATITUDE                                 *
C   RA_DAY , DEC_DAYDATE = CURRENT EPOCH RA AND DEC OF FIELD CENTRE    *
C   DAZ , DEL = AZ AND EL OFFSETS RELATIVE TO FIELD CENTRE      *
C
C Input:
C   DRA,DDEC   : Offset to the field centre (arcsec)
C   LST        : Local sideral time (sec)
C   SITE_LATITUDE     : Telescope latitude (degrees)
C   RA_DAY,DEC_DAY : RA and DEC of the day (degrees)
C Output:
C   DAZ,DEL   : AZ and EL offsets relative to the field centre (ARCSEC)
C******************************************************************
*
      SUBROUTINE NIC_RADEC_TO_AZEL(DRA,DDEC,LST,SITE_LATITUDE,RA_DAY,
     $DEC_DAY,DAZ,DEL)
*
      INCLUDE 'const.inc'
*
      REAL*8   DRA,DDEC
      REAL*8   LST
      REAL*8   SITE_LATITUDE
      REAL*8   RA_DAY,DEC_DAY
      REAL*8   DAZ,DEL
      REAL*8   HA
      REAL*8   SIN_DEC, COS_DEC
      REAL*8   SIN_DRA,COS_DRA,SIN_DDEC,COS_DDEC
      REAL*8   SIN_EL, COS_EL, SIN_PA,  COS_PA
      REAL*8   SIN_DAZ,COS_DAZ, SIN_DEL,  COS_DEL
      REAL*8   SIN_LAT,COS_LAT
      REAL*8   SIN_HA ,COS_HA
      REAL*8   PA,EL
*
      HA = -(LST/240.D0 - RA_DAY)*DEG_TO_RAD   ! hour angle
*
      SIN_LAT = DSIN(SITE_LATITUDE*DEG_TO_RAD)
      COS_LAT = DCOS(SITE_LATITUDE*DEG_TO_RAD)
      SIN_DEC = DSIN(DEC_DAY*DEG_TO_RAD)
      COS_DEC = DCOS(DEC_DAY*DEG_TO_RAD)
      SIN_HA  = DSIN(HA)
      COS_HA  = DCOS(HA)
* Compute Az, El corresponding to RA, DEC, and LST
      SIN_EL  = SIN_DEC*SIN_LAT + COS_DEC*COS_HA*COS_LAT
      EL      = DASIN(SIN_EL)
      COS_EL  = SQRT(1.D0-SIN_EL**2)
      SIN_PA  = COS_LAT*SIN_HA/COS_EL
      COS_PA  = (SIN_LAT - SIN_DEC*SIN_EL)/(COS_EL*COS_DEC)
      PA      = ATAN2(SIN_PA,COS_PA)
* Compute Az and El offsets (Daz and Del)
      SIN_DRA  = DSIN(DRA*SEC_TO_RAD)
      COS_DRA  = DCOS(DRA*SEC_TO_RAD)
      SIN_DDEC = DSIN(DDEC*SEC_TO_RAD)
      COS_DDEC = DCOS(DDEC*SEC_TO_RAD)
      SIN_DEL = -SIN_DRA*SIN_PA*COS_DDEC + COS_PA*SIN_DDEC
      DEL = DASIN(SIN_DEL)               ! del in ]-pi/2,pi/2]  (del~0)
      COS_DEL = DCOS(DEL)
      COS_DAZ = COS_DRA*COS_DDEC/COS_DEL
      IF (ABS(PA).GT.1.D-6) THEN
         SIN_DAZ = (-SIN_DDEC + COS_PA*SIN_DEL)/(SIN_PA*COS_DEL)
         DAZ = DATAN2(SIN_DAZ,COS_DAZ)
      ELSE
         DAZ = -DRA*SEC_TO_RAD
      ENDIF
* Conversion in arcsec
      DAZ = DAZ/SEC_TO_RAD
      DEL = DEL/SEC_TO_RAD
      RETURN
      END
C******************************************************************
*
c     NIC_AZEL_TO_RADEC
c
C   THIS ROUTINE TAKES AS INPUT THE AZ AND EL OFFSET RELATIVE     *
C   TO THE FIELD CENTRE, THE L.S.T., THE OBSERVATORY LATITUDE     *
C   AND THE CURRENT EPOCH EQUATORIAL COORDINATES OF THE FIELD     *
C   CENTRE AND RETURNS THE RA AND DEC OFFSETS RELATIVE TO THE     *
C   FIELD CENTRE.                                                 *
C   ALONG , ELAT = OFFSETS RELATIVE TO FIELD CENTRE IN AZ-EL      *
C   ALST = L.S.T. OF DATA POINT                                   *
C   TELLAT = OBSERVATORY LATITUDE                                 *
C   RADATE , DEDATE = CURRENT EPOCH RA AND DEC OF FIELD CENTRE    *
C   DRA , DDEC = RA AND DEC OFFSETS RELATIVE TO FIELD CENTRE      *
C
C Input:
C   ALONG,ALAT : Offset to the field centre (arcsec)
C   ALST       : Local sideral time (sec)
C   TELLAT     : Telescope latitude (degrees)
C   RADATE,DEDATE : RA and DEC of the day (degrees)
C Output:
C   DRA,DDEC   : RA and DEC offsets relative to the field centre (arcsec)
C******************************************************************
*
      SUBROUTINE NIC_AZEL_TO_RADEC(ALONG,ELAT,ALST,TELLAT,RADATE,
     $DEDATE,DRA,DDEC)
*
c      INCLUDE 'radec.inc'
      INCLUDE 'const.inc'
*
      REAL*8   ALONG
      REAL*8   ALST,ALST1
      REAL*8   ELAT,TELLAT
      REAL*8   RADATE,DEDATE,RACUR
      REAL*8   DRA,DDEC
      REAL*8   DLO,DLA,HA
      REAL*8   SINR,COSR,SINE,COSE,SINPA,COSPA,SINLO,COSLO,SINLA,COSLA
      REAL*8   SINDD,COSDD,COSDRA,SINDRA
      REAL*8   PA
      REAL*8   DECCUR,OBSLAT,SINDEC,COSDEC,SINOL,COSOL
*
*
      ALST1 = ALST / 240. * DEG_TO_RAD
      RACUR = RADATE * DEG_TO_RAD
      DECCUR = DEDATE * DEG_TO_RAD
      OBSLAT = TELLAT * DEG_TO_RAD
*
      DLO = ALONG / 3600. * DEG_TO_RAD
      DLA = ELAT /3600. * DEG_TO_RAD
      HA = RACUR - ALST1                 ! hour angle
*
      SINDEC = DSIN(DECCUR)              ! computed in RDCON
      COSDEC = DSQRT(1.D0 - SINDEC * SINDEC)
      SINOL = DSIN(OBSLAT)
      COSOL = DSQRT(1.D0 - SINOL * SINOL)
*
      SINR = DSIN(HA)
      COSR = DCOS(HA)
      SINE = SINDEC * SINOL + COSDEC * COSOL * COSR
C       E = DASIN(SINE)
      COSE = DSQRT(1.D0 - SINE * SINE)
      SINPA = COSOL * SINR / COSE
      COSPA = (SINOL - SINDEC *SINE) / (COSDEC * COSE)
      PA = DATAN2(SINPA,COSPA)
      SINLO = DSIN(DLO)
      COSLO = DCOS(DLO)
      SINLA = DSIN(DLA)
      COSLA = DCOS(DLA)
      SINDD = COSPA * SINLA - SINPA * COSLA *SINLO
      DDEC = DASIN(SINDD)
      COSDD = DCOS(DDEC)
      COSDRA = COSLA * COSLO / COSDD
      IF(ABS(SINPA).LT.1.D-6) THEN
         DRA = -DLO
      ELSE
         SINDRA = -(SINLA - COSPA * SINDD) / (SINPA * COSDD)
         DRA = DATAN2(SINDRA,COSDRA)
      ENDIF
      DRA = DRA / SEC_TO_RAD             ! DRA  in arcsecs
      DDEC = DDEC / SEC_TO_RAD           ! DDEC in arcsecs
      RETURN
      END
*
C ----------------------------------------------------------------------
C
C     SAVE_NMB()
C
C ----------------------------------------------------------------------
      SUBROUTINE SAVE_NMB(NMB_NAME,SUBSCAN_EPOCH,SUBSCAN_TIME,
     $LST,AZ,EL,XCOOR,YCOOR,SIG)
      INCLUDE 'nmb.inc'
      INCLUDE 'nic.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'par.inc'
      INCLUDE 'const.inc'
* parameters
      CHARACTER*(*)NMB_NAME
      REAL         SUBSCAN_EPOCH(1), SUBSCAN_TIME(1)
      REAL         LST(1),AZ(1),EL(1),XCOOR(1),YCOOR(1),SIG(1)
* local parameters
      INTEGER*4 IP_OUTPUTSTREAM,POINTER,IP,IADDR
      INTEGER*4 OUTPUTSTREAM
      INTEGER      LENC,NFILE
      INTEGER      I
      INTEGER      NUMBER, ILAST, LEN_RECORD
      INTEGER      ICHAN
      INTEGER      SIC_GETVM, IER
      CHARACTER*80    FILENAME
* close nmb file
      IF (INPUTSTREAM.NE.0) THEN
         CALL CLOSE_STREAM(INPUTSTREAM)
         INPUTSTREAM = 0
      ENDIF
*
* open file
      NFILE = LENC(NMB_NAME)
      FILENAME = NMB_NAME
      FILENAME(NFILE+1:NFILE+1) = CHAR(0)
      CALL OPEN_NEW_STREAM(FILENAME,OUTPUTSTREAM)
      IP_OUTPUTSTREAM=POINTER(OUTPUTSTREAM,MEMORY)
* save header
      HEADER_PREHEADER(5) =
     $'>FORMAT LST AZ EL SCAN_COORD SIGNAL'//CHAR(13)//CHAR(0)
      HEADER_PREHEADER(7)=
     $'>CHANNELS    '//CHAR(13)//CHAR(0)
      WRITE(HEADER_PREHEADER(7)(11:13),'(i3)') NCHAN
      HEADER_CHANNELS = NCHAN
      CALL WRITEHEADERINFORMATION(HEADER_STRUCT,MEMORY(IP_OUTPUTSTREAM))
* subscans
      ILAST = 0
      DO I=1,NSUBSCAN
         NUMBER = I
         IF (I.EQ.NSUBSCAN) ILAST = 1
         CALL WRITESUB(NUMBER,SUBSCAN_EPOCH(I),SUBSCAN_TIME(I),
     $   ILAST, MEMORY(IP_OUTPUTSTREAM))
      ENDDO
* records
*
* Record section
* ---------------
      LEN_RECORD = 5 + NCHAN             ! lst,az,el,2*scan_coord + n channels
      IF (INDEX(HEADER_PREHEADER(5),'<').NE.0) THEN
         LEN_RECORD = 5+2*NCHAN
      ENDIF
      IER = SIC_GETVM(LEN_RECORD,IADDR)
      IP = POINTER(IADDR,MEMORY)
*
      ILAST = 0
      DO I=1,NRECORD
         CALL R4TOR4(LST(I),MEMORY(IP),1)
         CALL R4TOR4(AZ(I), MEMORY(IP+1),1)
         CALL R4TOR4(EL(I), MEMORY(IP+2),1)
         CALL R4TOR4(XCOOR(I),MEMORY(IP+3),1)
         CALL R4TOR4(YCOOR(I),MEMORY(IP+4),1)
         DO ICHAN=1,NCHAN
            CALL R4TOR4(SIG(I+(ICHAN-1)*NRECORD),
     $      MEMORY(IP+5+ICHAN-1),1)
         ENDDO
         IF (I.EQ.NRECORD) ILAST = 1
         CALL WRITERECORD(MEMORY(IP),NCHAN,ILAST,
     $   MEMORY(IP_OUTPUTSTREAM))
      ENDDO
*
      CALL FREE_VM(LEN_RECORD,IADDR)
      CALL CLOSE_STREAM(OUTPUTSTREAM)
      END
