* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* onoff_flux.f
*
*-----------------------------------------------------------------------*
C+HPUX
c$OPTIMIZE LEVEL1
C-HPUX
      SUBROUTINE  ONOFF_FLUX (SIGNAL,NEW_SIGNAL,SIGWEI,FLUX_ONF,
     $FLUX_RMS,
     $INDX,NSUBREC,MEAN_ON,MEAN_OFF,X0,Y0,Z0,X1,Y1,
     $YSUBS,YSUBS_SIG,FLAG_SUBSCAN,SNOISEM,SNOISEW,CHAN_FLAG,OA,OE,
     $X,Y,Z,W,SEL_CHAN,PHASES,FITTING_METHOD,SIG_FIT,RMS_FIT,DIST,
     $SOURCE_CHAN,XCOO,YCOO,EL,NFIT)
*
      INCLUDE 'const.inc'
      INCLUDE 'onoff.inc'
      INCLUDE 'parameter.inc'
*
      REAL      SIGNAL(1),NEW_SIGNAL(1),SIGWEI(1)
      REAL      FLUX_ONF(NCHAN),FLUX_RMS(NCHAN)
      INTEGER   INDX(NSUBSCAN), NSUBREC(NSUBSCAN)
      REAL      MEAN_ON(NCHAN),MEAN_OFF(NCHAN)
      REAL      X0(NPAIRS),Y0(NPAIRS),Z0(NPAIRS)
      REAL      X1(NCHAN),Y1(NCHAN)
      REAL      YSUBS(NSUBSCAN,NCHAN),YSUBS_SIG(NSUBSCAN,NCHAN)
      LOGICAL   FLAG_SUBSCAN(NSUBSCAN,NCHAN)
      REAL      SNOISEM(NRECORD,3),SNOISEW(NRECORD)
      LOGICAL   CHAN_FLAG(NCHAN)
      REAL      OA(1),OE(1)
      REAL      X(1),Y(1),Z(1),W(1)
      INTEGER   SEL_CHAN(1)
      REAL      PHASES(NRECORD,2)
      INTEGER   FITTING_METHOD(NCHAN)
      INTEGER   NFIT
      REAL      SIG_FIT(NRECORD),RMS_FIT(NFIT),DIST(NCHAN,NCHAN)
      LOGICAL   SOURCE_CHAN(NCHAN)
      INTEGER   LIST(100),N
      REAL      XCOO(NRECORD),YCOO(NRECORD),EL(NRECORD)
* local variables
      INTEGER   ISUBSCAN,ICHAN,JCHAN
      REAL      A,B,C
      INTEGER   I,II,NP,NPH,IX
      INTEGER   NPH_ON, NPH_OFF
      REAL      WYMON,WYMOFF,MOY
      REAL      DIST_MAX
      INTEGER   IND
      REAL      DUMMY
      INTEGER   J,K,IFIT,BEST_FIT
      REAL      F(5)
      REAL      RMS_ON,RMS_OFF,RMS_MIN
      REAL      COS_EL,SIN_EL
      REAL      XWOB(2),XX,YY,HWOB
      REAL      MAXSIZE
*
* Initialize constants
      MAXSIZE = ONOFF_BEAM_RADIUS**2
      DIST_MAX = FIT_DIST
*
* Compute distances
      DO ICHAN=1,NCHAN
         DO JCHAN=1,NCHAN
            DIST(ICHAN,JCHAN) = SQRT((OA(ICHAN)-OA(JCHAN))**2
     $      +(OE(ICHAN)-OE(JCHAN))**2)
         ENDDO
      ENDDO
*
* Check if source is seen by channel
      XWOB(1) = XCOO(1)/10.
      XWOB(2) = (XCOO(1+NSUBREC(1)))/10.
      HWOB   = ABS(WOBBLER_THROW/2.)
      DO  ICHAN=1,NCHAN
         SOURCE_CHAN(ICHAN) = .FALSE.    ! default initialization
      ENDDO
      DO  ICHAN=1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            IF (ICHAN.EQ.REF_CHAN) THEN
               SOURCE_CHAN(ICHAN) = .TRUE.
            ELSE
               SOURCE_CHAN(ICHAN) = .FALSE.
               COS_EL = COS(FIELD_ROTATION*EL(NRECORD/2)/180.*PI)
               SIN_EL = SIN(FIELD_ROTATION*EL(NRECORD/2)/180.*PI)
               DO I=1,2
                  XX= XWOB(I)
     $            +(OA(REF_CHAN)-OA(ICHAN))*COS_EL
     $            +(OE(REF_CHAN)-OE(ICHAN))*SIN_EL
                  YY=
     $            -(OA(REF_CHAN)-OA(ICHAN))*SIN_EL
     $            +(OE(REF_CHAN)-OE(ICHAN))*COS_EL
                  IF (((XX-HWOB)**2+YY**2).LT.MAXSIZE)
     $            SOURCE_CHAN(ICHAN) = .TRUE.
                  IF (((XX+HWOB)**2+YY**2).LT.MAXSIZE)
     $            SOURCE_CHAN(ICHAN) = .TRUE.
               ENDDO
            ENDIF
            IF ( SOURCE_CHAN(ICHAN)) WRITE(6,*) 'Channel ' ,ICHAN,
     $      ' is in the beam defined by ONOFF_BEAM variable'
         ENDIF
      ENDDO
*
* Determine the begining index of each subscan
      INDX(1) = 1
      DO  ISUBSCAN = 2,NSUBSCAN
         INDX(ISUBSCAN) = INDX(ISUBSCAN-1) + NSUBREC(ISUBSCAN-1)
      ENDDO
*
* Calculate the mean of the signal subscan by subscan
      DO I =1,NRECORD
         SIGWEI(I)=1
      ENDDO
*
      IF (COMPUTE.EQ.1) THEN
*                                     ! compute weighted mean
         DO ICHAN =1,NCHAN
            DO ISUBSCAN=1,NSUBSCAN
               IX = INDX(ISUBSCAN)
               CALL WMEAN(SIGNAL(IX+(ICHAN-1)*NRECORD),
     $         SIGWEI,NSUBREC(ISUBSCAN),
     $         YSUBS(ISUBSCAN,ICHAN),YSUBS_SIG(ISUBSCAN,ICHAN),
     $         DUMMY)
            ENDDO
         ENDDO
      ELSE
         DO ICHAN =1,NCHAN
            DO ISUBSCAN=1,NSUBSCAN
               IX = INDX(ISUBSCAN)
               CALL MEDIAN(SIGNAL(IX+(ICHAN-1)*NRECORD),
     $         SIGWEI,NSUBREC(ISUBSCAN),
     $         YSUBS(ISUBSCAN,ICHAN),YSUBS_SIG(ISUBSCAN,ICHAN))
            ENDDO
         ENDDO
         DO I =1,NRECORD                 ! restore sigwei
            SIGWEI(I)=1
         ENDDO
      ENDIF
*
* Set chan_flag  if YSUBS(ISUBSCAN,ICHAN) equal to 0 for all subscans
      DO ICHAN =1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            DO ISUBSCAN=1,NSUBSCAN
               IF (YSUBS(ISUBSCAN,ICHAN).NE.0) GO TO 10
            ENDDO
10          IF (ISUBSCAN.GT.NSUBSCAN) THEN
               WRITE(6,*) ' *** No signal for channel ',ICHAN, ' *** '
               CHAN_FLAG(ICHAN)=.TRUE.
            ENDIF
         ENDIF
      ENDDO
*
* Initialize flag_subscan(nsubscan,nchan)
*
      DO ICHAN =1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            N = 0
            DO ISUBSCAN=1,NSUBSCAN
               IF (YSUBS(ISUBSCAN,ICHAN).EQ.0) THEN
                  FLAG_SUBSCAN(ISUBSCAN,ICHAN) = .TRUE.
                  N = N+1
                  LIST(N) = ISUBSCAN
               ELSE
                  FLAG_SUBSCAN(ISUBSCAN,ICHAN) = .FALSE.
!                Test if every point of subscan is blanked
                  DO I=INDX(ISUBSCAN),
     $               INDX(ISUBSCAN)+NSUBREC(ISUBSCAN)-1
                     II = I+(ICHAN-1)*NRECORD
                     IF (ABS(SIGNAL(II)-BLANKING).GT.0.5) GOTO 100
                  ENDDO
100               IF (I.EQ.INDX(ISUBSCAN)+NSUBREC(ISUBSCAN))THEN
                     FLAG_SUBSCAN(ISUBSCAN,ICHAN) = .TRUE.
                     N = N+1
                     LIST(N) = ISUBSCAN
                  ENDIF
               ENDIF
            ENDDO
            IF (N.NE.0) THEN
               WRITE(6,*)
     $         'I-ONOFF, Channel ',ICHAN,' flagged subscan(s): ',
     $         (LIST(I),I=1,N)
            ENDIF
         ENDIF
      ENDDO
*
* For unsymetrical scans, calculate and remove a first order baseline on
* signal and y_subs
*
      IF (ONOFF_MODE.EQ.UNSYMMETRICAL) THEN    ! only for unsymmetrical scans
         DO ICHAN =1,NCHAN
            IF (.NOT.CHAN_FLAG(ICHAN)) THEN
               DO NP=1,NPAIRS
                  NPH_ON = 2*NP-1
                  NPH_OFF = 2*NP
                  X0(NP) = 2*NP-0.5
                  Y0(NP) = (YSUBS(NPH_ON,ICHAN)+YSUBS(NPH_OFF,ICHAN))/2
                  Z0(NP) = 0
                  IF (.NOT.FLAG_SUBSCAN(NPH_ON,ICHAN) .AND.
     $            .NOT.FLAG_SUBSCAN(NPH_OFF,ICHAN)) THEN
                     Z0(NP) = 1/YSUBS_SIG(NPH_ON,ICHAN)**2+
     $               1/YSUBS_SIG(NPH_OFF,ICHAN)**2
                  ENDIF
               ENDDO
               CALL FITLINE(X0,Y0,Z0,NPAIRS, A, B)
               DO ISUBSCAN = 1,NSUBSCAN
                  YSUBS(ISUBSCAN,ICHAN) = YSUBS(ISUBSCAN,ICHAN)
     $            - (A*ISUBSCAN+B)
                  DO I = INDX(ISUBSCAN),
     $               INDX(ISUBSCAN)+NSUBREC(ISUBSCAN)-1
                     II = I+(ICHAN-1)*NRECORD
                     IF (ABS(SIGNAL(II)-BLANKING).GE.0.5) THEN
                        SIGNAL(II) = SIGNAL(II)-(A*ISUBSCAN+B)
                     ENDIF
                  ENDDO
               ENDDO
            ENDIF
         ENDDO
      ENDIF
*
*
* Subtraction of noise
* --------------------
      IF (NOISE) THEN
*
* Calculate the mean of the positive and negative phases
         DO ICHAN = 1,NCHAN
            IF (.NOT.CHAN_FLAG(ICHAN)) THEN
               MEAN_ON(ICHAN)= 0
               WYMON  = 0
               MEAN_OFF(ICHAN) = 0
               WYMOFF =  0
               DO NP =1,NPAIRS
                  IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN! ON-OFF symmetrical
                     IF (MOD(NP,2).EQ.1) THEN
                        NPH_ON  = 2*NP-1
                        NPH_OFF = 2*NP
                     ELSE
                        NPH_ON  = 2*NP
                        NPH_OFF = 2*NP-1
                     ENDIF
                  ELSE                   ! ON-OFF unsymmetrical
                     NPH_ON  = 2*NP-1
                     NPH_OFF = 2*NP
                  ENDIF
                  IF (.NOT.FLAG_SUBSCAN(NPH_ON,ICHAN)) THEN
                     MEAN_ON(ICHAN)= MEAN_ON(ICHAN) +
     $               YSUBS(NPH_ON,ICHAN)/YSUBS_SIG(NPH_ON,ICHAN)**2
                     WYMON  = WYMON + 1/YSUBS_SIG(NPH_ON,ICHAN)**2
                  ENDIF
                  IF (.NOT.FLAG_SUBSCAN(NPH_OFF,ICHAN)) THEN
                     MEAN_OFF(ICHAN) = MEAN_OFF(ICHAN) +
     $               YSUBS(NPH_OFF,ICHAN)/YSUBS_SIG(NPH_OFF,ICHAN)**2
                     WYMOFF = WYMOFF + 1/YSUBS_SIG(NPH_OFF,ICHAN)**2
                  ENDIF
               ENDDO
               IF (WYMON.NE.0) THEN
                  MEAN_ON(ICHAN)  = MEAN_ON(ICHAN) /WYMON
               ENDIF
               IF (WYMOFF.NE.0) THEN
                  MEAN_OFF(ICHAN) = MEAN_OFF(ICHAN)/WYMOFF
               ENDIF
            ENDIF
         ENDDO
*
* High frequency skynoise reduction by substracting the mean of adjacent
* channels
*
         DO ICHAN = 1,NCHAN
            IF (.NOT.CHAN_FLAG(ICHAN)) THEN
               DO IFIT=1,NFIT
                  DO I=1,NRECORD
                     SNOISEM(I,IFIT) = 0
                     SNOISEW(I) = 0
                  ENDDO
               ENDDO
*
               IF (FIT_TYPE.EQ.1.OR.FIT_TYPE.EQ.4) THEN
                  IFIT = 1
                  DO JCHAN=1,NCHAN
                     IF (JCHAN.NE.ICHAN .AND. JCHAN.NE.REF_CHAN) THEN
                        IF (.NOT.CHAN_FLAG(JCHAN)) THEN
                           IF (DIST(ICHAN,JCHAN).LE.DIST_MAX.AND.
     $                     .NOT.SOURCE_CHAN(JCHAN)) THEN
                              DO NPH= 1,NSUBSCAN
                                 IF (.NOT.FLAG_SUBSCAN(NPH,JCHAN)) THEN
                                    IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN
                                       IF (MOD(NPH,4).LE.1) THEN
                                          MOY = MEAN_ON(JCHAN)
                                       ELSE
                                          MOY = MEAN_OFF(JCHAN)
                                       ENDIF
                                    ELSE ! ON-OFF unsymmetrical
                                       IF (MOD(NPH,2).EQ.1) THEN
                                          MOY = MEAN_ON(JCHAN)
                                       ELSE
                                          MOY = MEAN_OFF(JCHAN)
                                       ENDIF
                                    ENDIF
                                    DO I=INDX(NPH),
     $                                 INDX(NPH)+NSUBREC(NPH)-1
                                       II = I+(JCHAN-1)*NRECORD
                                       IF (ABS(SIGNAL(II)-BLANKING)
     $                                 .GE.0.5) THEN
                                          SNOISEM(I,IFIT)=
     $                                    SNOISEM(I,IFIT)
     $                                    +(SIGNAL(II) - MOY)
     $                                    /YSUBS_SIG(NPH,JCHAN)**2
                                          SNOISEW(I)=SNOISEW(I)+
     $                                    1/YSUBS_SIG(NPH,JCHAN)**2
                                       ENDIF
                                    ENDDO
                                 ENDIF   ! ysubs_sig(nph,jchan).ne.0
                              ENDDO
                           ENDIF         ! if (dist.lt.distmax)
                        ENDIF
                     ENDIF
                  ENDDO
*
                  DO I=1,NRECORD
                     IF (SNOISEW(I).GT.0) THEN
                        SNOISEM(I,IFIT) = SNOISEM(I,IFIT)/SNOISEW(I)
                     ENDIF
                  ENDDO
               ENDIF
*
* High frequency skynoise reduction by computing plane around considered
* channel
*
               IF (FIT_TYPE.EQ.2.OR.FIT_TYPE.EQ.4) THEN
                  IFIT = 2
*     Find 6 channels around considered channel
                  IND = 1
                  DO JCHAN=1,NCHAN
                     IF (JCHAN.NE.REF_CHAN.AND..NOT.CHAN_FLAG(JCHAN)
     $               .AND.JCHAN.NE.ICHAN) THEN
                        IF (DIST(ICHAN,JCHAN).LE.DIST_MAX.AND.
     $                  .NOT.SOURCE_CHAN(JCHAN)) THEN
                           SEL_CHAN(IND) = JCHAN
                           X(IND) = OA(JCHAN)
                           Y(IND) = OE(JCHAN)
                           IND = IND+1
                        ENDIF
                     ENDIF
                  ENDDO
                  IND = IND-1
c                 write(6,*)ichan,':',(sel_chan(i),i=1,ind)
                  IF (IND.GE.3) THEN
                     DO ISUBSCAN=1,NSUBSCAN
                        DO I=1,IND
                           W(I)=0
                           IF (.NOT.FLAG_SUBSCAN(ISUBSCAN,SEL_CHAN(I)))
     $                     W(I) = 1/YSUBS(ISUBSCAN,SEL_CHAN(I))**2
                        ENDDO
                        DO J=INDX(ISUBSCAN),
     $                     INDX(ISUBSCAN)+NSUBREC(ISUBSCAN)-1
                           DO I=1,IND
                              MOY = MEAN_ON(SEL_CHAN(I))
                              IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN
                                 IF (MOD(ISUBSCAN,4).EQ.2.OR.
     $                           MOD(ISUBSCAN,4).EQ.3)
     $                           MOY = MEAN_OFF(SEL_CHAN(I))
                              ELSE
                                 IF (MOD(ISUBSCAN,2).EQ.0)
     $                           MOY = MEAN_OFF(SEL_CHAN(I))
                              ENDIF
                              II = J+(SEL_CHAN(I)-1)*NRECORD
                              Z(I) = SIGNAL(II) - MOY
                           ENDDO
                           CALL LINEAR_FIT_PLANE(X,Y,Z,W,IND,A,B,C)
                           SNOISEM(J,IFIT)=
     $                     A*OA(ICHAN)+B*OE(ICHAN)+C
                        ENDDO
                     ENDDO
                  ENDIF
               ENDIF
*
* High frequency skynoise reduction by computing paraboloid around considered
* channel
*
               IF (FIT_TYPE.EQ.3.OR.FIT_TYPE.EQ.4) THEN
                  IFIT=3
* Find n channels around reference channel
                  IND = 1
                  DO JCHAN=1,NCHAN
                     IF (JCHAN.NE.REF_CHAN.AND..NOT.CHAN_FLAG(JCHAN)
     $               .AND.JCHAN.NE.ICHAN) THEN
                        IF (DIST(ICHAN,JCHAN).LE.DIST_MAX.AND.
     $                  .NOT.SOURCE_CHAN(JCHAN)) THEN
                           SEL_CHAN(IND) = JCHAN
                           X(IND) = OA(JCHAN)
                           Y(IND) = OE(JCHAN)
                           IND = IND+1
                        ENDIF
                     ENDIF
                  ENDDO
                  IND = IND-1
c                  write(6,*)ichan,ind,':',(sel_chan(i),i=1,ind)
                  IF (IND.GE.5) THEN
                     DO ISUBSCAN=1,NSUBSCAN
                        DO I=1,IND
                           W(I) = 0
                           IF (.NOT.FLAG_SUBSCAN(ISUBSCAN,SEL_CHAN(I)))
     $                     W(I) = 1/YSUBS(ISUBSCAN,SEL_CHAN(I))**2
                        ENDDO
                        DO J=INDX(ISUBSCAN),INDX(ISUBSCAN)
     $                     +NSUBREC(ISUBSCAN)-1
                           DO I=1,IND
                              MOY = MEAN_ON(SEL_CHAN(I))
                              IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN
                                 IF (MOD(ISUBSCAN,4).EQ.2.OR.
     $                           MOD(ISUBSCAN,4).EQ.3)
     $                           MOY = MEAN_OFF(SEL_CHAN(I))
                              ELSE
                                 IF (MOD(ISUBSCAN,2).EQ.0)
     $                           MOY = MEAN_OFF(SEL_CHAN(I))
                              ENDIF
                              II = J+(SEL_CHAN(I)-1)*NRECORD
                              Z(I) = SIGNAL(II) - MOY
                           ENDDO
                           CALL FIT_QUADRATIC(X,Y,Z,W,IND,F)
                           SNOISEM(J,IFIT)=
     $                     F(1)+F(2)*OA(ICHAN)+F(3)*OE(ICHAN)+
     $                     F(4)*OA(ICHAN)*OA(ICHAN)+
     $                     F(5)*OE(ICHAN)*OE(ICHAN)
                        ENDDO
                     ENDDO
                  ENDIF
               ENDIF
*
* If fit_method "auto", compute rms to find the best fitting method
               IF (FIT_TYPE.EQ.4) THEN
                  DO IFIT=1,NFIT
                     DO I=1,NRECORD
                        II = I +(ICHAN-1)*NRECORD
                        SIG_FIT(I)=SIGNAL(II)-SNOISEM(I,IFIT)
                     ENDDO
                     J = 1
                     K = 1
                     DO NP= 1,NPAIRS
                        IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN
                           IF (MOD(NP,2).EQ.1) THEN
                              NPH_ON  = 2*NP-1
                              NPH_OFF = 2*NP
                           ELSE
                              NPH_ON  = 2*NP
                              NPH_OFF = 2*NP-1
                           ENDIF
                        ELSE             ! ON-OFF unsymmetrical
                           NPH_ON  = 2*NP-1
                           NPH_OFF = 2*NP
                        ENDIF
                        IF (.NOT.FLAG_SUBSCAN(NPH_ON,ICHAN).AND.
     $                  .NOT.FLAG_SUBSCAN(NPH_OFF,ICHAN)) THEN
                           DO I=INDX(NPH_ON),INDX(NPH_ON)+
     $                        NSUBREC(NPH_ON)-1
                              IF (ABS(SIG_FIT(I)-BLANKING).GE..5) THEN
                                 PHASES(J,1)=SIG_FIT(I)
                                 J = J+1
                              ENDIF
                           ENDDO
                           DO I=INDX(NPH_OFF),INDX(NPH_OFF)
     $                        +NSUBREC(NPH_OFF)-1
                              IF (ABS(SIG_FIT(I)-BLANKING).GE.5) THEN
                                 PHASES(K,2)=SIG_FIT(I)
                                 K = K+1
                              ENDIF
                           ENDDO
                        ENDIF
                     ENDDO
                     J = J-1
                     K = K-1
                     CALL WMEAN(PHASES(1,1),SIGWEI,J,DUMMY,DUMMY,RMS_ON)
                     CALL WMEAN(PHASES(1,2),SIGWEI,K,DUMMY,DUMMY,
     $               RMS_OFF)
                     RMS_FIT(IFIT) = (RMS_ON+RMS_OFF)/2/SQRT(2.)
                  ENDDO
                  RMS_MIN=1.E32
                  DO IFIT=1,NFIT
                     IF (RMS_FIT(IFIT).LT.RMS_MIN) THEN
                        RMS_MIN = RMS_FIT(IFIT)
                        BEST_FIT= IFIT
                     ENDIF
                  ENDDO
               ENDIF
*
* Recalculate the mean of the signal after skynoise correction
* for each subscan by subscan.
               FITTING_METHOD(ICHAN) = FIT_TYPE
               IF (FIT_TYPE.EQ.4) FITTING_METHOD(ICHAN)  = BEST_FIT
               DO I= 1,NRECORD
                  II = I+(ICHAN-1)*NRECORD
                  IF (ABS(SIGNAL(II)-BLANKING).GE.0.5) THEN
                     NEW_SIGNAL(II) = SIGNAL(II)
     $               - SNOISEM(I,FITTING_METHOD(ICHAN))
                  ELSE
                     NEW_SIGNAL(II) = BLANKING
                  ENDIF
               ENDDO
*
            ENDIF                        ! endif (.not.chan_flag)
         ENDDO                           ! do ichan=1,nchan
*
* Copy NEW_SIGNAL in SIGNAL at the end of noise correction calculation
         DO ICHAN=1,NCHAN
            IF (.NOT.CHAN_FLAG(ICHAN)) THEN
               DO I=1,NRECORD
                  II = I+(ICHAN-1)*NRECORD
                  SIGNAL(II) = NEW_SIGNAL(II)
               ENDDO
               DO ISUBSCAN =1,NSUBSCAN
                  IX = INDX(ISUBSCAN)
                  CALL WMEAN(SIGNAL(IX+(ICHAN-1)*NRECORD),
     $            SIGWEI,NSUBREC(ISUBSCAN),
     $            YSUBS(ISUBSCAN,ICHAN),
     $            YSUBS_SIG(ISUBSCAN,ICHAN),DUMMY)
               ENDDO
            ENDIF                        ! endif (.not.chan_flag)
         ENDDO
      ENDIF                              ! endif (noise)
*
*
* Calculate the flux (with or without noise suppression)
*
      DO ICHAN=1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            DO NP= 1,NPAIRS
               IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN   ! ON-OFF symmetrical
                  IF (MOD(NP,2).EQ.1) THEN
                     NPH_ON  = 2*NP-1
                     NPH_OFF = 2*NP
                  ELSE
                     NPH_ON  = 2*NP
                     NPH_OFF = 2*NP-1
                  ENDIF
               ELSE                      ! ON-OFF unsymmetrical
                  NPH_ON  = 2*NP-1
                  NPH_OFF = 2*NP
               ENDIF
               X0(NP) = (YSUBS(NPH_ON,ICHAN)-YSUBS(NPH_OFF,ICHAN))/2.  ! signal
               IF (.NOT.FLAG_SUBSCAN(NPH_ON,ICHAN) .AND.
     $         .NOT.FLAG_SUBSCAN(NPH_OFF,ICHAN)) THEN
                  Y0(NP) = 1/YSUBS_SIG(NPH_ON,ICHAN)**2+
     $            1/YSUBS_SIG(NPH_OFF,ICHAN)**2! weight
               ELSE
                  Y0(NP) = 0
               ENDIF
            ENDDO
            CALL WMEAN(X0, Y0, NPAIRS, FLUX_ONF(ICHAN),
     $      DUMMY,FLUX_RMS(ICHAN))
         ELSE
            FLUX_ONF(ICHAN) = 0
            FLUX_RMS(ICHAN) = 0
         ENDIF
      ENDDO
*
* Compute global RMS
      DO ICHAN=1,NCHAN
         J = 1
         K = 1
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            DO NP= 1,NPAIRS
               IF (ONOFF_MODE.EQ.SYMMETRICAL) THEN   ! ON-OFF symmetrical
                  IF (MOD(NP,2).EQ.1) THEN
                     NPH_ON  = 2*NP-1
                     NPH_OFF = 2*NP
                  ELSE
                     NPH_ON  = 2*NP
                     NPH_OFF = 2*NP-1
                  ENDIF
               ELSE                      ! ON-OFF unsymmetrical
                  NPH_ON  = 2*NP-1
                  NPH_OFF = 2*NP
               ENDIF
               IF (.NOT.FLAG_SUBSCAN(NPH_ON,ICHAN).AND.
     $         .NOT.FLAG_SUBSCAN(NPH_OFF,ICHAN)) THEN
                  DO I=INDX(NPH_ON),INDX(NPH_ON)+NSUBREC(NPH_ON)-1
                     II = I+(ICHAN-1)*NRECORD
                     IF (ABS(SIGNAL(II)-BLANKING).GE..5) THEN
                        PHASES(J,1)=SIGNAL(II)
                        J = J+1
                     ENDIF
                  ENDDO
                  DO I=INDX(NPH_OFF),INDX(NPH_OFF)+
     $               NSUBREC(NPH_OFF)-1
                     II = I+(ICHAN-1)*NRECORD
                     IF (ABS(SIGNAL(II)-BLANKING).GE..5) THEN
                        PHASES(K,2)=SIGNAL(II)
                        K = K+1
                     ENDIF
                  ENDDO
               ENDIF
            ENDDO
            J = J-1
            K = K-1
            CALL WMEAN(PHASES(1,1),SIGWEI,J,DUMMY,DUMMY,RMS_ON)
            CALL WMEAN(PHASES(1,2),SIGWEI,K,DUMMY,DUMMY,RMS_OFF)
            FLUX_RMS(ICHAN) = (RMS_ON+RMS_OFF)/2/SQRT(2.)
         ENDIF
      ENDDO
*
*  Compute weighted mean on adjacent channels
*
      IND = 1
      NFLG = 0
      X1(1) = 0
      DO ICHAN=1,NCHAN
         IF (ICHAN.NE.REF_CHAN.AND..NOT.CHAN_FLAG(ICHAN)) THEN
            IF (DIST(ICHAN,REF_CHAN).LE.DIST_MAX.AND.
     $      .NOT.SOURCE_CHAN(ICHAN)) THEN
               IF (.NOT.CHAN_FLAG(ICHAN)) THEN
                  X1(IND) = FLUX_ONF(ICHAN)
                  Y1(IND) = 1/FLUX_RMS(ICHAN)**2
                  IND = IND+1
               ELSE
                  NFLG = NFLG+1
               ENDIF
            ENDIF
         ENDIF
      ENDDO
      IND = IND-1
      IF (IND.GT.1) THEN
         CALL WMEAN(X1, Y1, IND, ADJACENT_FLUX,DUMMY,RMS_ADJACENT_FLUX)
      ELSE
         ADJACENT_FLUX = X1(1)
         RMS_ADJACENT_FLUX =0
      ENDIF
      REF_FLUX = FLUX_ONF(REF_CHAN)-ADJACENT_FLUX
      REF_FLUX_RMS = FLUX_RMS(REF_CHAN)**2+RMS_ADJACENT_FLUX**2
      REF_FLUX_RMS = SQRT(REF_FLUX_RMS)
      NADJ = IND
      END
*******************************************************************
*
* weighted mean
* signal(npts) is the input variable to mean
* weight(npts) is the input variable for the weights
* npts is the number of data to consider
* mean is the estimated mean
* sigma  is the standard deviation
* sigmam is the estimated error
*
      SUBROUTINE WMEAN(SIGNAL,WEIGHT,NPTS,MEAN,SIGMA,SIGMAM)
*
      INCLUDE 'parameter.inc'            ! blanking value
      INTEGER  NPTS,NP
      REAL     SIGNAL(NPTS),WEIGHT(NPTS)
      REAL     MEAN,SIGMA,SIGMAM
*
      INTEGER  I
      REAL*8     XMEAN,XSIGMA,XSIGMAM,WEISUM
*
* Calculate the mean
      IF (NPTS.LE.1) THEN
         WRITE(6,*) 'Any missing points ...'
         MEAN = SIGNAL(1)
         SIGMA  = 0
         SIGMAM = 0
      ELSE
         XMEAN = 0
         WEISUM = 0
         XSIGMA = 0
         NP = 0
         DO I = 1,NPTS
            IF (ABS(SIGNAL(I)-BLANKING).GE.0.5) THEN
               XMEAN = XMEAN + WEIGHT(I)*SIGNAL(I)
               WEISUM = WEISUM + WEIGHT(I)
               NP = NP+1
            ENDIF
         ENDDO
         IF (WEISUM.NE.0) THEN
            XMEAN = XMEAN / WEISUM
         ELSE
c            WRITE(6,*) 'No data available (blanked ?) !'
            XMEAN= 0
         ENDIF
         DO I= 1,NPTS
            IF (ABS(SIGNAL(I)-BLANKING).GE.0.5) THEN
               XSIGMA = XSIGMA + WEIGHT(I)*(SIGNAL(I)-XMEAN)**2
            ENDIF
         ENDDO
         IF (WEISUM.NE.0.AND.NP.GT.1) THEN
            XSIGMA = DSQRT(XSIGMA*NP/((NP-1)*WEISUM))
            XSIGMAM = XSIGMA/DSQRT(DBLE(NP))
         ELSE
            XSIGMA = 0
            XSIGMAM=0
         ENDIF
         MEAN = XMEAN
         SIGMA= XSIGMA
         SIGMAM = XSIGMAM
      ENDIF
      END
*******************************************************************
*
* median mean
* &1 is the input variable to mean, its dimension is &2
* &2 is a work area, its dimension is &2
* &3 is the number of data to consider on the variable &1 and &2
* &4 is the estimated mean
* &5 is the standard deviation
*
*
      SUBROUTINE MEDIAN(SIGNAL,WORK,NPTS_IN,MEAN,DXMED)
*
      INCLUDE 'parameter.inc'
      INTEGER  NPTS_IN
      REAL     SIGNAL(NPTS_IN),WORK(NPTS_IN)
      REAL     MEAN,DXMED
*
      INTEGER  I,L,J,IR,N2,NPTS
      REAL     W
* Sort data Heap sort algorythm (Numerical Recipes p231)
      NPTS = 0
      DO I=1,NPTS_IN
         IF (ABS(SIGNAL(I)-BLANKING).GE.0.5) THEN
            WORK(I)=SIGNAL(I)            ! copy signal in work area
            NPTS = NPTS+1
         ENDIF
      ENDDO
      IF (NPTS.EQ.0) THEN
         WRITE(6,*)' *** MEDIAN: no valid data ***'
         MEAN = 0
         DXMED = 1
         RETURN
      ENDIF
*
      L = NPTS/2 + 1
      IR = NPTS
      DO WHILE (.TRUE.)
         IF (L.GT.1) THEN
            L = L-1
            W = WORK(L)
         ELSE
            W = WORK(IR)
            WORK(IR) = WORK(1)
            IR = IR-1
            IF (IR.EQ.1) THEN
               WORK(1) = W
               GOTO 100                  ! end
            ENDIF
         ENDIF
         I = L
         J = L+1
         DO WHILE (J.LE.IR)
            IF (J.LT.IR) THEN
               IF(WORK(J).LT.WORK(J+1)) J = J+1
            ENDIF
            IF (W.LT.WORK(J)) THEN
               WORK(I)=WORK(J)
               I = J
               J = J+J
            ELSE
               J = IR+1
            ENDIF
         ENDDO
         WORK(I) = W
      ENDDO
*
* Compute median mean
100   N2 = NPTS/2
      IF (2*N2 .EQ. NPTS) THEN
         MEAN = 0.5*(WORK(N2)+WORK(N2+1))
      ELSE
         MEAN = WORK(N2+1)
      ENDIF
*
* deviation
      DXMED = WORK( NINT (REAL(N2)+SQRT(REAL(NPTS))/2.) ) -
     $WORK( NINT (REAL(N2)-SQRT(REAL(NPTS))/2.) )
      DXMED = DXMED*0.5*SQRT(REAL(NPTS))
      IF (DXMED.EQ.0) DXMED = 1.
      END
*************************************************************************
* Procedure fitline.nic
* Calculate the coeff (A,B) of the baseline (Ax+B)
* &1 the variable x[npairs]
* &2 the variable y[npairs]
* &3 the respectives weights w[npairs]
* &4 the number of data = npairs
* &5 A
* &6 B
*
      SUBROUTINE FITLINE(X,Y,W,NPTS,A,B)
*
      INTEGER NPTS
      REAL    X(NPTS),Y(NPTS),W(NPTS),A,B
*
      REAL    D, SW, SWX, SWY, SWXY, SWX2
      INTEGER I
      IF (NPTS.GT.1) THEN
         SW = 0
         SWX = 0
         SWY = 0
         SWXY = 0
         SWX2 = 0
         DO I=1,NPTS
            SW  = SW + W(I)
            SWX = SWX+ W(I)*X(I)
            SWY = SWY+ W(I)*Y(I)
            SWXY = SWXY + W(I)*X(I)*Y(I)
            SWX2 = SWX2 + W(I)*X(I)**2
         ENDDO
         D = SW*SWX2-SWX**2
         IF (ABS(D).GT.1.D-32) THEN
            A = (SW*SWXY-SWX*SWY)/D
            B = (SWY*SWX2-SWX*SWXY)/D
         ELSE
            A = 0
            B = 0
         ENDIF
      ELSE
         A = 0
         B = Y(1)
      ENDIF
      END
*************************************************************************
* Procedure ONOFF_FLUX_SHARC
*
      SUBROUTINE  ONOFF_FLUX_SHARC (SIGNAL,FLUX_ONF,FLUX_RMS,
     $X0,Y0,
     $CHAN_FLAG)
*
      INCLUDE 'onoff.inc'
      INCLUDE 'parameter.inc'
*
      REAL      SIGNAL(1)
      REAL      FLUX_ONF(NCHAN),FLUX_RMS(NCHAN)
      REAL      X0(NPAIRS),Y0(NPAIRS)
      LOGICAL   CHAN_FLAG(NCHAN)
      INTEGER   ICHAN
      INTEGER   NP
      INTEGER   NPH_ON, NPH_OFF
      INTEGER   IOFF_CHAN
      REAL      DUMMY
*
      IF (NOISE) THEN
         WRITE(6,'(/,a,/)')
     $   '*** Noise reduction not available with SHARC data  ***'
         NOISE = .FALSE.
      ENDIF
*
*
*
* Calculate the flux
*
      DO ICHAN=1,NCHAN
         IF (.NOT.CHAN_FLAG(ICHAN)) THEN
            IOFF_CHAN=(ICHAN-1)*NRECORD
            DO NP= 1,NPAIRS
* ON-OFF symmetrical  (ON OFF OFF ON)
               IF (MOD(NP,2).EQ.1) THEN
                  NPH_ON  = 2*NP-1 +IOFF_CHAN
                  NPH_OFF = 2*NP   +IOFF_CHAN
               ELSE
                  NPH_ON  = 2*NP   +IOFF_CHAN
                  NPH_OFF = 2*NP-1 +IOFF_CHAN
               ENDIF
               IF (ABS(SIGNAL(NPH_ON)-BLANKING).GT.0.1
     $         .AND.ABS(SIGNAL(NPH_OFF)-BLANKING).GT.0.1) THEN
                  X0(NP) = (SIGNAL(NPH_ON)-SIGNAL(NPH_OFF))/2.
                  Y0(NP) = 1             ! weight
               ELSE
                  X0(NP) = BLANKING
                  Y0(NP) = 0
               ENDIF
            ENDDO
            CALL WMEAN(X0, Y0, NPAIRS, FLUX_ONF(ICHAN),
     $      DUMMY,FLUX_RMS(ICHAN))
         ELSE
            FLUX_ONF(ICHAN) = 0
            FLUX_RMS(ICHAN) = 0
         ENDIF
      ENDDO
      END
*************************************************************************
      SUBROUTINE LINEAR_FIT_PLANE (X,Y,Z,W,N,A,B,C)
!
! where z = a*x+b*y+c
!
      INTEGER I,N
      REAL*8 G,GX,GY,GZ,GXX,GYY,GXY,GZX,GZY
      REAL X(N),Y(N),Z(N),W(N)
      REAL A,B,C,D
*
      G   = 0
      GX  = 0
      GY  = 0
      GZ  = 0
      GXX = 0
      GYY = 0
      GXY = 0
      GZX = 0
      GZY = 0
      DO I=1,N
         G   = G   + W(I)
         GX  = GX  + W(I)*X(I)
         GY  = GY  + W(I)*Y(I)
         GZ  = GZ  + W(I)*Z(I)
         GXX = GXX + W(I)*X(I)*X(I)
         GYY = GYY + W(I)*Y(I)*Y(I)
         GXY = GXY + W(I)*X(I)*Y(I)
         GZX = GZX + W(I)*Z(I)*X(I)
         GZY = GZY + W(I)*Z(I)*Y(I)
      ENDDO
      D = G*(GXX*GYY-GXY*GXY)+GX*(GY*GXY-GX*GYY)+GY*(GX*GXY-GY*GXX)
      IF (D.NE.0) THEN
         A = GZ*(GY*GXY-GX*GYY)+GZX*(G*GYY-GY*GY)+GZY*(GX*GY-G*GXY)
         B = GZ*(GX*GXY-GY*GXX)+GZX*(GX*GY-G*GXY)+GZY*(G*GXX-GX*GX)
         C = GZ*(GXX*GYY-GXY*GXY)+GZX*(GY*GXY-GX*GYY)+
     $   GZY*(GX*GXY-GY*GXX)
         A = A/D
         B = B/D
         C = C/D
      ENDIF
      RETURN
      END
