********************************************************
*   RESULTP DE TSUKUTTA FILE WO YONDE FOLDING WO SURU
********************************************************
      REAL*8 STPD,STP,DTP
      DIMENSION D1(500),D2(500),R1(500),DAY(500),ELV(500)
      DIMENSION ED1(500),ED2(500),ER1(500)
      CHARACTER CMNT*128,FNAME*20,CMNT0(30)*128
      DIMENSION CHI(1000),FCURVE(200),FSIGM(200),NFC(200)
      logical exflag
      dimension ntime0(7),ntime1(7)
      real*8 mjdhjm,mjdowr
ce    character fname1*30
      character fname1*60

      ntmax = 1000
      DISCL=0
      DISCH=99999999.
      STTD=0
      STPD=99999999.
*-------- READ OUTPUT PARAMETERS -----------------------
*-------- read source file name ------------------------
 1000 continue
      write(*,*) 'Please input source file name !'
      read(*,'(a)') fname
      write(*,*)
ce--- write(fname1,'(''~egoshi/ASM.database/Source.data/'',a)') fname
      write(fname1,'(''/home/kitamoto/Asm/Asm/data/'',a)') fname
      inquire(file=fname1,exist=exflag)
      if(.not.exflag)then
*         write(fname1,'(''~egoshi/ASM.database/Source.data/'',a)')
*     &        fname
*         inquire(file=fname1,exist=exflag)
*         if (.not.exflag) then
         write(*,*) fname1
         write(*,*) 'File cannot be found ! Input once again!'
         write(*,*)
         goto 1000
*       endif
      endif
*-------- read start and stop day ---------------------
 1100 continue
ce---
 1201 do 1200 i12=4,7
         ntime0(i12)=0
         ntime1(i12)=0
 1200 continue
      write(*,*) 'Input Start day.'
      write(*,*) 
      write(*,*) 'Year or Day from 1987 Jan 1 ?'
      read(*,*) sday
      if (sday.lt.1987 .or. sday.gt.1992) then
         mjdhjm=dble(sday)+46796.0d0
         write(*,*)
         goto 1250
      endif
      ntime0(1)=sday
      write(*,*) 'Month ?'
      read(*,'(i2)') ntime0(2)
      write(*,*) 'Day ?'
      read(*,'(i2)') ntime0(3)
      write(*,*) 
      call mjulia(ntime0,mjdhjm)
 1250 continue
      write(*,*) 'Input Stop day.'
      write(*,*) 
      write(*,*) 'Year or Day from 1987 Jan 1 ?'
      read(*,*) sday
      if (sday.lt.1987 .or. sday.gt.1992) then
         mjdowr=dble(sday)+46796.0d0
         write(*,*)
         goto 1300
      endif
      ntime1(1)=sday
      write(*,*) 'Month ?'
      read(*,'(i2)') ntime1(2)
      write(*,*) 'Day ?'
      read(*,'(i2)') ntime1(3)
      write(*,*) 
      call mjulia(ntime1,mjdowr)
 1300 continue
ce---
      sttd=mjdhjm-46796.0d0
      stpd=mjdowr-46796.0d0
      IF(STPD.LE.0)STPD=99999999.
*-------- read energy range ----------------------------
      IHR=1
      write(*,*) 'Input energy range (channel)'
      write(*,*)
      write(*,*) 'Lower channel?'
      read(*,*) ich11
      write(*,*)
      write(*,*) 'Higher channel?'
      read(*,*) ich12
      write(*,*)
      ich1=ich11*100+ich12

      write(*,*) 'Input higher energy range for hardness ratio'
      write(*,*) '(input 0 ... donot plot hardness ratio)'
      write(*,*)
      write(*,*) 'Lower channel?'
      read(*,*) ich21
      write(*,*)
      IF(ICH21.EQ.0)then
         IHR=0
         goto 1400
      endif
      write(*,*) 'Higher channel?'
      read(*,*) ich22
      write(*,*)
      ich2=ich21*100+ich22
*-------- read disclimination level -------------------
 1400 continue
      write(*,*) 'Input lower disclimination level of flux'
      read(*,*) discl
      write(*,*)
      write(*,*) 'Input upper desclimination level of flux'
      read(*,*) dischd
      write(*,*)
      IF(DISCHD.GT.0)DISCH=DISCHD
*-------- read effective area -------------------------
ce---
*      write(*,*) 'Input relative effective area of detector'
*      read(*,*) area
*      write(*,*)
*      IF(AREA.LE.0.0)AREA=1.0
      area=1.0
ce---
*-------- read trial parameters for folding  ----------
      write(*,*) 'Input trial parameters for folding'
      write(*,*)
      write(*,*) 'Start period (day)'
      read(*,*) stp
      write(*,*)
      write(*,*) 'Incliment period (day)'
      read(*,*) dtp
 1600 continue
      write(*,*)
      write(*,*) 'Trial number of incliment'
      read(*,*) ntry
      if(ntry.gt.ntmax)then
         write(*,*) 'Trial number should be less than',ntmax
         go to 1600
      endif
      write(*,*)
      write(*,*) 'Bin number for plotting folding result'
      read(*,*) nbin
      write(*,*)
      IF(NBIN.GT.200)NBIN=200
*-------- FILE OPEN -----------------------------------
      OPEN(7,FILE=FNAME1,STATUS='old')
*-------- READ DATA -----------------------------------
      CALL READD
     &     (ICH1,ICH2,STTD,STPD,D1,ED1,D2,ED2,DAY,ELV,NDATA,
     &     CMNT,AREA,DISCL,DISCH)

*      print *,(d1(j),j=1,10)
      CLOSE (7)
ce
      CMNT=fname(1:13)//' '//CMNT
*------- SET HARDNESS RATIO ----------------------------
      IF(IHR.EQ.1)
     &     CALL HRATIO(D1,ED1,D2,ED2,NDATA,R1,ER1)
*------- FOLDING   -------------------------------------
      CALL FOLD
     &     (STP,DTP,NTRY,NBIN,CHI,FCURVE,FSIGM,NFC,D1,ED1,DAY
     &     ,NDATA,PER0)
*------- OUTPUT    -------------------------------------
      CALL FOLDO
     &     (STTD,STPD,STP,DTP,
     &     NTRY,NBIN,CHI,FCURVE,FSIGM,NFC,PER0,DISCL,DISCH,
     &     CMNT0,NCMNT)
*------- PLOT DATA -------------------------------------
* 9000 continue
*      write(*,*)
*      write(*,*) '------------------------'
*      write(*,*) '  1. Light curve'
*      write(*,*) '  2. Folding'
*      write(*,*) '  3. Try another period'
*      write(*,*)
*      write(*,*) '  7. Try again'
*      write(*,*) '  8. Another source'
*      write(*,*) '  9. End'
*      write(*,*) '------------------------'
*      write(*,*)
*      write(*,*) 'Please input the number'
*      write(*,*)
*      read(*,*) num
*      write(*,*)
*      if (num.eq.1) then
*         CALL DPLOT(ICH1,ICH2,D1,ED1,D2,ED2,R1,ER1,DAY,ELV,NDATA
*     &        ,CMNT,IHR)
*      elseif (num.eq.2) then
*         CALL DPLOT2
*     &        (D1,ED1,DAY,NDATA,PER0,STP,DTP,FCURVE,FSIGM,NBIN,
*     &        CHI,NTRY,ICH1,CMNT,CMNT0,NCMNT)
*      elseif (num.eq.3) then
*         write(*,*) 'Input period (day)'
*         read(*,*) per
*         write(*,*)
*         call fold2(fcurve,fsigm,nfc,per,nbin,day,d1,ed1,ndata)
*         CALL DPLOT2
*     &        (D1,ED1,DAY,NDATA,PER,STP,DTP,FCURVE,FSIGM,NBIN,
*     &        CHI,NTRY,ICH1,CMNT,CMNT0,NCMNT)
*      elseif (num.eq.7) then
*         goto 1100
*      elseif (num.eq.8) then
*         goto 1000
*      elseif (num.eq.9) then
*         goto 999
*      endif
*      goto 9000
*------- END -------------------------------------------
 999  CONTINUE
      STOP
      END
********************************************************
      SUBROUTINE READD
     &(ICH1,ICH2,STTD,STPD,D1,ED1,D2,ED2,DAY,ELV,NDATA,CMNT,
     &           AREA,DISCL,DISCH)
*-------------------------------------------------------
*     ASM NO DATA WO YOMU
*     INPUT PARAMETERS
*       ICH1,ICH2 : ENERGY CHANNEL
*       STTD,STPD : START AND STOP DAY
*       AREA      : MENSEKI
*       DISCL     : LOWER DISCLIMINATION LEVEL
*       DISCH     : UPPER DISCLIMINATION LEVEL
*     OUTPUT PARAMETERS
*       D1,ED1    : DATA1
*       D2,ED2    : DATA2
*       DAY       : DAY
*       ELV       : ELEVATION
*       NDATA     : DATA NO KAZU
*       CMNT      : CMMENT
*-------------------------------------------------------
      parameter (nmax=500)
      DIMENSION D1(nmax),D2(nmax),DAY(nmax),ELV(nmax)
      DIMENSION ED1(nmax),ED2(nmax)
      CHARACTER DUMMY*128
      CHARACTER CMNT*128
      REAL*8 STPD,STP,DTP
      integer nchk(nmax)
      real*4 data(nmax,19),edata(nmax,19)
      real*4 data1(19),edata1(19)
*------- clear data ------------
 2001 do 2000 i21=1,nmax
 2101    do 2100 i22=1,19
            data(i21,i22)=0.0
            edata(i21,i22)=0.0
 2100    continue
 2000 continue
*-------- READ COMMENT -------------------
      READ(7,'(A128)')CMNT
      NDATA=0
 1000 CONTINUE
*-------- READ DATA ----------------------
      READ(7,'(A128)',END=999)DUMMY
      READ(DUMMY,'(I4)')NO
      IF(NO.EQ.0)GO TO 1000
      READ(DUMMY,'(4X,I4,6F8.0)')ICH,D,E1,E2,ED,DA,EL

*      print *," read;",ich,d,e1,e2,ed,da,el
      IF(DA.LT.STTD.OR.DA.GT.STPD) GO TO 1000
csk
*      if(ich.eq.ich1 .or. ich.eq.ich2) then
*         IF(D.LE.DISCL) GO TO 1000
*         IF(D.GT.DISCH) GO TO 1000
*        IF(DA.LT.STTD.OR.DA.GT.STPD) GO TO 1000
*      endif
csk
      if(ich.eq.108) then
         mch=17
      elseif(ich.eq.916) then
         mch=18
      elseif(ich.eq.116) then
         mch=19
      else
         mch=int(float(ich)/101)
      endif
      IF(NDATA.LE.0)GO TO 1200
*------- SEARCH PREVIOUS DAY ------------
 1101 DO 1100 I11=1,NDATA
         IF(DAY(I11).EQ.DA)THEN
            data(i11,mch)=d
            edata(i11,mch)=ed
            elv(i11)=el
            GO TO 1000
         ENDIF
 1100 CONTINUE
*------- NEW DATA -----------------------
 1200 CONTINUE
      NDATA=NDATA+1
      data(ndata,mch)=d
      edata(ndata,mch)=ed
      day(ndata)=da
      ELV(NDATA)=EL
*------- CHECK NDATA---------------------
      IF(NDATA.GE.NMAX)THEN
         WRITE(6,'('' TOO MANY DATA '')')
         GO TO 999
      ENDIF
*------- NEXT DATA-----------------------
      GO TO 1000
 999  CONTINUE

*      print *,"data(1,j)",(data(1,j),j=1,19)
*      print *,"data(2,j)",(data(2,j),j=1,19)
*      print *,"data(3,j)",(data(3,j),j=1,19)
*      print *,"data(4,j)",(data(4,j),j=1,19)


*------- for all channel -------
      if(ndata.le.0) goto 990
 1301 do 1300 i13=1,ndata
         if(data(i13,19).eq.0.0) then
            if(data(i13,17).ne.0.0.and.data(i13,18).ne.0.0) then
               data(i13,19)=data(i13,17)+data(i13,18)
               edata(i13,19)=sqrt(edata(i13,17)**2
     &              +edata(i13,18)**2)
            else
               data(i13,19)=0.0
               edata(i13,19)=0.0
            endif
         endif
 1300 continue
ce------ Check bad data --------
 1401 do 1400 i14=1,ndata
         nchk(i14)=1
 1411    do 1410 i141=1,18
            if(data(i14,i141).le.-0.3) then
               nchk(i14)=0
               goto 1400
            elseif(edata(i14,i141).ge.0.5) then
               nchk(i14)=0
               goto 1400
            elseif(edata(i14,i141).le.2e-3) then
               nchk(i14)=0
               goto 1400
            elseif(data(i14,i141).le.0) then
               zero = data(i14,i141) + edata(i14,i141)*3
               if(zero.lt.0) then
                  nchk(i14)=0
                  goto 1400
               endif
            endif
 1410    continue
 1400 continue
*------- Sort data -------------
      mdata=0
 1531 do 1530 i153=1,ndata
         if(nchk(i153).eq.0) goto 1530
         mdata=mdata+1
         day(mdata)=day(i153)
 1541    do 1540 i154=1,19
            data(mdata,i154) =data(i153,i154)
            edata(mdata,i154)=edata(i153,i154)
 1540    continue
 1530 continue
      ndata=mdata
*------- gain correction ----------------
* 3001 do 3000 i3=1,ndata
*         print *,"before data(i3,1)=",data(i3,1)
* 3101    do 3100 i31=1,19
*            data1(i31)=data(i3,i31)
*            edata1(i31)=edata(i3,i31)
* 3100    continue
*         call crceng(day(i3),data1,edata1)
* 3201    do 3200 i32=1,19
*            data(i3,i32)=data1(i32)
*            edata(i3,i32)=edata1(i32)
* 3200    continue
*         print *,"after data(i3,1)=",data(i3,1)
* 3000 continue
*------- ich1 and ich2 data -------------
      if(ich1.eq.108) then
         mch1=17
         mch11=17
         mch12=17
      elseif(ich1.eq.916) then
         mch1=18
         mch11=18
         mch12=18
      elseif(ich1.eq.116) then
         mch1=19
         mch11=19
         mch12=19
      else
         mch11=int(float(ich1)/100)
         mch12=ich1-mch11*100
      endif
      if(ich2.eq.108) then
         mch2=17
         mch21=17
         mch22=17
      elseif(ich2.eq.916) then
         mch2=18
         mch21=18
         mch22=18
      elseif(ich2.eq.116) then
         mch2=19
         mch21=19
         mch22=19
      else
         mch21=int(float(ich2)/100)
         mch22=ich2-mch21*100
      endif
 4001 do 4000 i4=1, ndata
         d1(i4)=0
         ed1(i4)=0
*         print *,i4," data 1= ",d1(i4),ed1(i4),  ich1, ich2, mch11,mch12
 4101    do 4100 i41=mch11,mch12
            d1(i4)=d1(i4)+data(i4,i41)/area
            ed1(i4)=ed1(i4)+(edata(i4,i41)/area)**2
 4100    continue
         ed1(i4)=sqrt(ed1(i4))

*         print *,i4," data 1= ",d1(i4),ed1(i4)

         if (ich2.eq.0) goto 4000
         d2(i4)=0
         ed2(i4)=0
 4201    do 4200 i42=mch21,mch22
            d2(i4)=d2(i4)+data(i4,i42)/area
            ed2(i4)=ed2(i4)+(edata(i4,i42)/area)**2
 4200    continue
         ed2(i4)=sqrt(ed2(i4))
*         print *,i4,"data 2=",d2(i4),ed2(i4)
 4000 continue

*------- Sort data -------------
      mdata=0
 5001 do 5000 i5=1,ndata
         IF(d1(i5).LE.DISCL) GO TO 5000
         IF(d1(i5).GT.DISCH) GO TO 5000
         mdata=mdata+1
         day(mdata)=day(i5)
         d1(mdata) =d1(i5)
         ed1(mdata)=ed1(i5)
         d2(mdata) =d2(i5)
         ed2(mdata)=ed2(i5)
 5000 continue

*      print *,"in readd",(d1(j),j=1,10)
      ndata=mdata

      print *,"The number of data points are ",ndata

*------- RETURN -------------------------
 990  continue
      RETURN
      END
********************************************************
      SUBROUTINE HRATIO(D1,ED1,D2,ED2,NDATA,R1,ER1)
*-------------------------------------------------------
*     DATA NO RATIO WO KEISAN SURU
*     INPUT AREGUMENTS
*           D1,ED1  : DATA1
*           D2,ED2  : DATA2
*           NDATA   : DATA NO KAZU
*     OUTPUT AREGUMENTS
*           R1(I) =D2(I)/D1(I)
*           ER1(I) =SQRT((D1(I)*ED2(I)+D2(I)*ED2(I))/D2(I)**2)
*-------------------------------------------------------
      DIMENSION D1(500),D2(500),R1(500)
      DIMENSION ED1(500),ED2(500),ER1(500)
 1001 DO 1000 I1=1,NDATA
        IF(D2(I1).LE.0.OR.D1(I1).LE.0)THEN
          R1(I1)=-999.
          ER1(I1)=-999.
        ELSE
          R1(I1) =D2(I1)/D1(I1)
          ER1(I1) =SQRT((D2(I1)**2*ED1(I1)**2+
     &                   D1(I1)**2*ED2(I1)**2))/D1(I1)**2
        ENDIF
 1000 CONTINUE
      RETURN
      END
********************************************************
      SUBROUTINE FOLD
     &  (STP,DTP,NTRY,NBIN,CHI,FCURVE,FSIGM,NFC,D1,ED1,DAY,NDATA,PER0)
      DIMENSION CHI(NTRY),FCURVE(NBIN),FSIGM(NBIN),NFC(NBIN)
      DIMENSION D1(NDATA),DAY(NDATA),ED1(NDATA)
*-------------------------------------------------------
*     FOLDING
*        STP  : START PERIOD (DAY)
*        DTP  : INCLIMENT PERIOD (DAY)
*        NTRY : THE NUMBER OF INCLIMENT
*        CHI  : THE ARRAY OF CHI SQUARE
*        FCURVE : LIGHT CURVE AT THE MAXIMUM CHISQUARED PERIOD
*        FSIGM  : HYOUJUNN HENSA
*        NFC    : DATA NO KAZU
*        D1   : DATA
*        DAY  : DAY
*-------------------------------------------------------
        REAL*8 STPD,STP,DTP
*------- SEARCH CHI SQUARED ----------------------
 1001   DO 1000 I1=1,NTRY
*------- SET PERIOD ------------------------------
          PER=STP+DTP*(I1-1)
          CALL FOLD2(FCURVE,FSIGM,NFC,PER,NBIN,DAY,D1,ED1,NDATA)
*------- SET CHISQUARED---------------------------
          CALL CHIS(FCURVE,FSIGM,NFC,NBIN,CHI(I1))
*          print *," check in fold",i1,per, chi(i1)
 1000   CONTINUE
*------- SET MAX CHISQUARED---------------------------
        MC=1
 2001   DO 2000 I2=2,NTRY
          IF(CHI(MC).LT.CHI(I2))MC=I2
 2000   CONTINUE
        PER=STP+DTP*(MC-1)
        CALL FOLD2(FCURVE,FSIGM,NFC,PER,NBIN,DAY,D1,ED1,NDATA)
        PER0=PER
        RETURN
        END
**************************************************
        SUBROUTINE FOLD2(FCURVE,FSIGM,NFC,PER,NBIN,DAY,D1,ED1,NDATA)
        DIMENSION FCURVE(NBIN),FSIGM(NBIN),NFC(NBIN)
csk
        dimension fcurv2(200)
csk
        DIMENSION DAY(NDATA),D1(NDATA),ED1(NDATA)
*------- SET DPER   ------------------------------
          DPER=PER/FLOAT(NBIN)
*          print *," dper=",dper," ndata=",ndata," in fold2"
*------- CLEAR DATA ------------------------------
 1101     DO 1100 I11=1,NBIN
            FCURVE(I11)=0
            FSIGM(I11) =0
csk
            fcurv2(i11)=0
csk
            NFC(I11)   =0
 1100     CONTINUE
*------- ADD DATA   ------------------------------
 1201     DO 1200 I12=1,NDATA
*             print *,i12, day(i12),d1(i12)
            if (ed1(i12).le.0) goto 1200
            SIGMA=1/ED1(I12)**2
ce          N12=MOD(DAY(I12)/DPER,NBIN)+1
ce          N12=MOD(int(DAY(I12)/DPER),NBIN)+1
            n12 = (day(i12)/per - int(day(i12)/per)) * nbin + 1
csk         FCURVE(N12)=FCURVE(N12)+D1(I12)*SIGMA
csk
            FCURVE(N12)=FCURVE(N12)+D1(I12)
            fcurv2(n12)=fcurv2(n12)+d1(i12)**2
csk
            FSIGM(N12) =FSIGM(N12)+SIGMA
            NFC(N12)   =NFC(N12)+1
 1200     CONTINUE
 1251     DO 1250 I125=1,NBIN
            IF(NFC(I125).EQ.0) GO TO 1250
ce
            if(nfc(i125).eq.1) then
               fsigm(i125)=1/sqrt(fsigm(i125))
               goto 1250
            endif
csk         FCURVE(I125)=FCURVE(I125)/FSIGM(I125)
csk         FSIGM(I125)=1/SQRT(FSIGM(I125))

            FCURVE(I125)=FCURVE(I125)/float(nfc(i125))
            fcurv2(i125)=fcurv2(i125)/float(nfc(i125))
            FSIGM(I125)=SQRT((fcurv2(i125)-fcurve(i125)**2)
     &           /float(nfc(i125)))


csk

 1250     CONTINUE

c***       debug *****
* 2001     do 2000 i2=1, nbin
*             print *,i2," ",fcurve(i2)," ",fsigm(i2)
* 2000     continue


          RETURN
          END
**************************************************
       SUBROUTINE CHIS(FCURVE,FSIGM,NFC,NBIN,CHI)
       DIMENSION FCURVE(NBIN),FSIGM(NBIN),NFC(NBIN)
*-------------------------------------------------
*        SET CHISQUARED
*-------------------------------------------------
       A=0
       B=0
ce
       nfc0=0

 1001  DO 1000 I1=1,NBIN
ce       IF(NFC(I1).LE.0)GO TO 1000
         IF(NFC(I1).LE.0) then
            nfc0=nfc0+1
            goto 1000
         endif
         A=A+(1./FSIGM(I1)**2)
         B=B+(FCURVE(I1)/FSIGM(I1)**2)
 1000  CONTINUE
       RMEAN=B/A
       C=0
 1101  DO 1100 I11=1,NBIN
         IF(NFC(I11).LE.0)GO TO 1100
         C=C+((FCURVE(I11)-RMEAN)/FSIGM(I11))**2
 1100  CONTINUE
ce     CHI=C
       CHI=C/(float(nbin)-nfc0-1.0)
       RETURN
       END
********************************************************
      SUBROUTINE FOLDO
     & (STTD,STPD,STP,DTP,
     &  NTRY,NBIN,CHI,FCURVE,FSIGM,NFC,PER0,DISCL,DISCH,CMNT0,NCMNT)
      DIMENSION CHI(NTRY),FCURVE(NBIN),FSIGM(NBIN),NFC(NBIN)
      CHARACTER CMNT0(30)*128
      REAL*8 STPD,STP,DTP
*-------------------------------------------------------
*      OUTPUT THE RESULTS OF FOLDING
*-------------------------------------------------------
      write(*,*)
      WRITE(6,'(''*** RESULTS OF FOLDING ****'' )')
      WRITE(6,'(''  * DATA START   :'',F10.3,''DAY'')')STTD
      WRITE(6,'(''  * DATA STOP    :'',F10.3,''DAY'')')STPD
      WRITE(6,'(''  * START PERIOD :'',F10.3,''DAY'')')STP
      WRITE(6,'(''  * INCLIMENT    :'',F10.4,''DAY'')')DTP
      WRITE(6,'(''  * TRY          :'',I5)')NTRY
      WRITE(6,'(''  * BIN          :'',I5)')NBIN
      WRITE(6,'(''  * RESULT PERIOD:'',E15.7,''DAY'')') PER0
      WRITE(6,'(''  * LOWER DISCL  :'',E15.7)') DISCL
      WRITE(6,'(''  * UPPER DISCL  :'',E15.7)') DISCH
      write(*,*)
      WRITE(CMNT0(1),'(''*** RESULTS OF FOLDING ****'' )')
      WRITE(CMNT0(2),'(''  * DATA START   :'',F10.3,''DAY'')')STTD
      WRITE(CMNT0(3),'(''  * DATA STOP    :'',F10.3,''DAY'')')STPD
      WRITE(CMNT0(4),'(''  * START PERIOD :'',F10.3,''DAY'')')STP
      WRITE(CMNT0(5),'(''  * INCLIMENT    :'',F10.4,''DAY'')')DTP
      WRITE(CMNT0(6),'(''  * TRY          :'',I5)')NTRY
      WRITE(CMNT0(7),'(''  * BIN          :'',I5)')NBIN
      WRITE(CMNT0(8),'(''  * RESULT PERIOD:'',E15.7,''DAY'')')
     &     PER0
      WRITE(CMNT0(9),'(''  * LOWER DISCLI :'',E15.7)')DISCL
      WRITE(CMNT0(10),'(''  * UPPER DISCLI :'',E15.7)')DISCH
      NCMNT=10
ce--- 
      WRITE(6,'(''*** CHI SQUARE SEARCH ***''/
     &     ,''  NO.  PERIOD(DAY)   CHISQUARED  '')')
 1001 DO 1000 I1=1,NTRY
        PER=STP+DTP*(I1-1)
        WRITE(6,'(I6,F15.6,E15.7)')I1,PER,CHI(I1)
 1000 CONTINUE
ce---
      WRITE(6,'(''*** FOLDED LIGHT CURVE  ***''/
     &     ''   NO.   #      AVERAGE      S.D         '')')
 2001 DO 2000 I2=1,NBIN
         WRITE(6,'(I6,I4,F15.6,E15.7)')
     &     I2,NFC(I2),FCURVE(I2),FSIGM(I2)
 2000 CONTINUE
      write(*,*)
      RETURN
      END
****************************************************************
        subroutine mjulia(TIME,MJD)
        INTEGER*4  TIME(7)
        REAL*8     MJD
C
C  DATE AND TIME ----> MJD
C
C     INPUT DATA; TIME(7)
C                  YEAR, MONTH, DAY, HOUR, MINUTE, SEC, MILISEC
C                  IF YEAR < 100  , YEAR IS ASSUME TO BE 1900'S
C     OUT PUT DATA; MJD
C                  MJD=JD-2400000.5
C            DAY OF WEEK = MOD(MJD+3,7)  (0=SUN, 1=MON, ETC)
C
        REAL*8 Y,M,D
        IF(TIME(2).GT.2) THEN
          Y=DFLOAT(TIME(1))
          M=DFLOAT(TIME(2))
         ELSE
          Y=DFLOAT(TIME(1)-1)
          M=DFLOAT(TIME(2)+12)
        ENDIF
        IF(TIME(1).LT.100) Y=Y+1900.0
        D=DFLOAT(TIME(3))+
     &    (DFLOAT(TIME(4))+(DFLOAT(TIME(5))+
     &    (DFLOAT(TIME(6))+DFLOAT(TIME(7))/1000.0)/60.0 )/60.0) 
     &	   /24.0
        MJD=DINT(Y*365.25)+DINT(Y/400.)-DINT(Y/100.)
     &       +DINT((M-2.)*30.59) + D -678912.0D0
        RETURN
        END
c=====================================================================
      subroutine crceng(day,data,edata)
c---------------------------------------------------------------------
*      implicit NONE
c  input
      real*4 day,data(19),edata(19)
c  var
ce--- real*4 engch1(16)/1.000,1.709,2.255,2.861,
*    &                  3.466,4.084,4.737,5.329,
*    &                  5.967,7.124,8.320,9.505,
*    &                  10.698,12.979,15.393,17.898/
*     real*4 engch2(16)/1.709,2.255,2.861,3.466,
*    &                  4.084,4.737,5.329,5.967,
*    &                  7.124,8.320,9.505,10.698,
ce---&                  12.979,15.393,17.898,20.270/

*      real*4 engch1(16)/1.000,1.640,2.152,2.717,
*     &                  3.281,3.857,4.469,5.024,
*     &                  5.634,6.743,7.899,9.055,
*     &                  10.223,12.468,14.860,17.369/
*      real*4 engch2(16)/1.640,2.152,2.717,3.281,
*     &                  3.857,4.469,5.024,5.634,
*     &                  6.743,7.899,9.055,10.223,
*     &                  12.468,14.860,17.369,19.827/

	real*4 engch1(16), engch2(16)
      real*4 data2(19),edata2(19),egcrc1(16),egcrc2(16)
      real*4 gain,engch0,factor,factor2
      integer i1,i11,i2,i21,i3,ich
c  common
      common /engrng/ engch1,engch2
c  function
      real*4 calgin
      real*4 effic
c  begin


 1001 do 1000 i1=1,19
         data2(i1)=0
         edata2(i1)=0
 1000 continue

c  calculate gain
      gain = calgin(day)

 1101 do 1100 i11=1,16
         egcrc1(i11) = engch1(i11) / gain
         egcrc2(i11) = engch2(i11) / gain
 1100 continue

c  correction of each enegy ch
 2001 do 2000 i2=1,16

         if ( i2 .eq. 1 ) then

            factor = effic(engch1(1),engch1(2))
     &             / effic(egcrc1(1),egcrc1(2))
            data2(i2) = factor * data(1)
            edata2(i2) = factor * edata(1)
            goto 2000
         elseif ( i2 .eq. 2 ) then
            factor = effic(engch1(2),egcrc2(1))
     &             / effic(egcrc1(1),egcrc2(1))
            factor2 = effic(egcrc1(2),engch2(2))
     &              / effic(egcrc1(2),egcrc2(2))
            data2(i2) = factor * data(1) + factor2 * data(2)
            edata2(i2) = sqrt( ( factor * edata(1) )**2
     &                       + ( factor2 * edata(2) )**2 )
            goto 2000
         endif

         engch0 = engch1(i2)
         ich = 0
 2101    do 2100 i21=1,16
            if ( engch1(i2) .ge. egcrc1(i21)
     &           .and. engch1(i2) .lt. egcrc2(i21) ) then
               ich = i21
               goto 2200
            endif
 2100    continue

 2200    continue
         if ( engch2(i2) .le. egcrc2(ich) ) then
            factor = ( engch2(i2) - engch0 ) 
     &               / ( egcrc2(ich) - egcrc1(ich) )
            data2(i2) = data2(i2) + factor * data(ich)
            edata2(i2) = sqrt( edata2(i2)**2
     &                      + ( factor*edata(ich) )**2 )
         else
            factor = ( egcrc2(ich) - engch0 )
     &               / ( egcrc2(ich) - egcrc1(ich) )
            data2(i2) = data2(i2) + factor * data(ich)
            edata2(i2) = sqrt( edata2(i2)**2
     &                      + ( factor*edata(ich) )**2 )
            engch0 = egcrc2(ich)
            ich = ich + 1
            goto 2200
         endif

 2000 continue

c  correction of enegy ch 1-8, 9-16, 1-16
      factor = ( engch2(8)-engch1(1) ) / ( egcrc2(8)-egcrc1(1) )
      data2(17) = factor * data(17)
      edata2(17) = factor * edata(17)

      factor = ( egcrc2(8)-engch1(9) ) / ( egcrc2(8)-egcrc1(1) )
      factor2 = ( engch2(16)-egcrc1(9) ) / ( egcrc2(16)-egcrc1(9) )
      data2(18) = factor * data(17) + factor2 * data(18)
      edata2(18) = sqrt( ( factor*edata(17) )**2
     &                + ( factor2*edata(18) )**2 )

      factor = ( engch2(16)-engch1(1) ) / ( egcrc2(16)-egcrc1(1) )
      data2(19) = factor * data(19)
      edata2(19) = factor * edata(19)

 3001 do 3000 i3=1,19
         data(i3) = data2(i3)
         edata(i3) = edata2(i3)
 3000 continue


      return
      end
c=====================================================================
      function calgin(day)
c---------------------------------------------------------------------
      implicit NONE
c  input
      real*4 day
c  output
      real*4 calgin
c  const
      real*4 P1, P2, P3, P4
      parameter ( P1=7.849, P2=-3.1320e-4,
     &            P3=-2.1046e-3, P4=-2.1426e-4 )
c  var
      real*4 peak, offset
c  begin
      if ( day .lt. 1300.0 ) then
         peak = P1 + P2*day
      elseif ( day .lt. 1400.0 ) then
         peak = P1 + P2*1300.0 + P3*(day-1300.0)
      else
         peak = P1 + P2*1300.0 + P3*100.0 + P4*(day-1400.0)
      endif
      offset = P1 / ( 5.896 - 1.000 )
      calgin = ( peak + offset ) / ( P1 + offset )
      return
      end
c=====================================================================
      function effic(e1,e2)
c---------------------------------------------------------------------
      implicit NONE
c  input
      real*4 e1,e2
c  output
      real*4 effic
c  function
      real*4 eff,e,eff1,eff5,p1,p2
      eff(e) = 1.0 - exp( -p1 * e**p2 )
      eff1 = 1.0 - 0.001
      eff5 = 1.0 - 0.9
      p1 = -log(eff1)
      p2 = log( -log(eff5) / p1 ) / log(5.0)
c  begin

      if ( e2 .lt. e1 ) then
         write(*,*) 'ERROR e2 < e1 in function effic.'
         stop
      endif
      effic = ( eff(e1) + eff(e2) ) * ( e2 - e1 ) / 2.0

      return
      end

***********************************************************
	subroutine setengch

	real*4 engch1(16)
	real*4 engch2(16)
      common /engrng/ engch1,engch2

        engch1(1)=1.000
        engch1(2)=1.640
        engch1(3)=2.152   
        engch1(4)=2.717
        engch1(5)=3.281
        engch1(6)=3.857
        engch1(7)=4.469
        engch1(8)=5.024
        engch1(9)=5.634
        engch1(10)=6.743
        engch1(11)=7.899
        engch1(12)=9.055
        engch1(13)=10.223
        engch1(14)=12.468
        engch1(15)=14.860
        engch1(16)=17.369

        engch2(1)=1.640
        engch2(2)=2.152
        engch2(3)=2.717
        engch2(4)=3.281
        engch2(5)=3.857
        engch2(6)=4.469
        engch2(7)=5.024
        engch2(8)=5.634
        engch2(9)=6.743
        engch2(10)=7.899
        engch2(11)=9.055
        engch2(12)=10.223
        engch2(13)=12.468
        engch2(14)=14.860
        engch2(15)=17.369
        engch2(16)=19.827

	return
	end
