!!!!!!!!!!!!!!!!!!!!!
!!  Change History !!
!!!!!!!!!!!!!!!!!!!!!
!   TODO:
!   - [ ] Add a test for the DM[1-2][R-W] stuff in the python when we're happy with it
!
! <<< Known "Bugs" >>>
!   - Some models in the  ID_TABLE  (a.k.a mentioned in the input files) are not in the other tables of the
!   database. Likely cause: They can't be opened for ones reason or another.
!
! v2.0.2 - 7 Apr 2022
!   - eve_maker renamed to eve
!   - Removed the DELTA_* columns from the fortran - now only recording DM[1-2][R-W] WITH CONVERSION TO MSOL/YR 
!   - Added avgDM[1-2][R-W]_msolpyr to allow comparison of RLOF Vs wind contirbution to mass loss 
!      -> Specific mass loss value can be done on a case by case basis in python. It requires thorough
!      -> consideration of the bin edges (in the middle of each age step) to calculate the dt. 
!      -> WARNING: Do not use the factor of 1.989 that is mentioned in the manual. 
!   - Case determination no longer includes ONe core as edge cases have been found 
!     with sporadic ONe cores prior to CO core formation
!
! v2.0.1 - 4 Apr 2022
!   -  README.md  updated
!   - Merged [PR#6](https://github.com/UoA-Stars-And-Supernovae/EvE/pull/6) - FIX Case A bug
!   - Added a script and two files in the tests folder to help test that. These are temporary files and will be deleted.
!
! v2.0.0 - 3 Apr 2022
!   Created eve-v2 branch. Contains:
!   - eve_maker.f08 which creates text files of each table. Current tables implemented in the Schema are:
!     ID_TABLE, SUMMARY, CEE, MT, DEATH
!   - A python script eve_txt_to_hdf.py which combines the csv tables into an hdf5 file
!   - A test folder with test_eve.py which checks that the models in  SUMMARY  table are in `ID_TABLE` and vice versa,
!   that all the models in SUMMARY are also in `DEATH` (and vice versa), that all the models flagged with CEE or MT in
!   SUMMARY  also are recorded in the `CEE` and `MT` tables (and vice versa).
!   - A  README.md  has been added
!


program eve
implicit none

integer :: DEBUG_MODE
integer :: i, j ! i=loop iterations, j = index for the dummy arrays
integer :: modeltype
real*8, dimension(59) :: dummy
real*8 :: dump
real*8 :: modelimf, mixedimf, mixedage, weight
real*8 :: initblackhole, initperiod
real*8 :: m1_ZAMS, m2_ZAMS, period_ZAMS

character(len=10) :: imf_name, zname
character(len=300) :: models_dirname, input_dirname, input_file
character(len=100) :: fname_sneplot
character(len=4) :: znames(13)

! FILE I/O
integer :: rlfo_rows, n, m, err_status, err_status_sneplot, err_status_output! if gets value other than 0 we'll know somethi
character(256) :: err_iomsg
logical filecheck

! TIMING THE CODE
real ::  t1,t2

character(len=3) :: sn_type

integer :: MODEL_ID, TIMESTEP_start, TIMESTEP_end, MT_bool, CEE_bool, MT_bool_summary, CEE_bool_summary
real*8 :: AGE_start, M1_start, M2_start, P_start, LOGA_start, lifetime_MS
real*8 :: AGE_end, M1_end, M2_end, P_end, LOGA_end
real*8 :: AGE_previous, avgDM1R_msolpyr, avgDM2R_msolpyr, avgDM1W_msolpyr, avgDM2W_msolpyr  
character(len=17) ::  case_str, start_rlof_case

logical :: He_core_bool, CO_core_bool

! For each BPASSS Metallicity
DEBUG_MODE = 2
! INPUT FILE READ
models_dirname='../BPASS_hoki_dev/bpass-v2.2-newmodels/'
input_dirname = '../BPASS_hoki_dev/bpass_v2.2.1_imf135_300/'
imf_name='imf135_300'
zname='z010'
znames=(/'zem5', 'zem4','z001', 'z002', 'z003', 'z004','z006', 'z008', 'z010','z014', 'z020', 'z030', 'z040'/)

MODEL_ID=-1

do i=1,13,1
    err_status=0
    zname=znames(i)
    write(*,*) zname
    input_file=trim(input_dirname)//"input_bpass_"//trim(zname)//'_bin_'//trim(imf_name)
    if(DEBUG_MODE.eq.2) write(*,*) 'READING THROUGH--- ', input_file

    open(200, FILE=input_file,&
              STATUS='old',&
              ACCESS='SEQUENTIAL',&
              ACTION='read')

    read(200,*) dump ! this should be 1 but used to be a useful value - we don't use it so i called it dump

    ! SNE PLOT FORMAT
    211     format(59E16.7)

    ! Open file that will contain the SUMMARY table
    open(100, file="./eve_txt/"//trim(zname)//"_ID_TABLE.dat",&
              status='unknown', &
              access='sequential',&
              action='write')

    open(110, file="./eve_txt/"//trim(zname)//"_SUMMARY.dat",&
               status='unknown', &
               access='sequential',&
               action='write')

    write(100, *) "MODEL_ID ", "filenames"
    ! ID TABLES (MODEL ID and filenames)
    101     format(i8,'  ', a85)

    write(110, *) "MODEL_ID ","M1_ZAMS ", "M2_ZAMS ", "P_ZAMS ", "lifetime_MS ","MT_bool ","CEE_bool ",&
                  "AGE_end ", "M1_end ", "M2_end ", "P_end ","modelimf ","mixedimf ","type"
    ! Format statement SUMMARY table
    111     format(i8, 4E16.7, 2i2, 6E16.7, i3)

    ! Open file that will contain the SUMMARY table
    open(120, file="./eve_txt/"//trim(zname)//"_CEE.dat",&
              status='unknown', &
              access='sequential',&
              action='write')

    open(122, file="./eve_txt/"//trim(zname)//"_MT.dat",&
              status='unknown', &
              access='sequential',&
              action='write')

    ! Format statement CEE and MT table
    write(120, *) 'MODEL_ID ', 'TIMESTEP_start ','AGE_start ','M1_start ','M2_start ','P_start ','LOGA_start ',&
                  'TIMESTEP_end ','AGE_end ','M1_end ','M2_end ','P_end ','LOGA_end ',&
                  'avgDM1R_msolpyr ','avgDM2R_msolpyr ',&
                  'avgDM1W_msolpyr ', 'avgDM2W_msolpyr ', 'CASE '

    write(122, *) 'MODEL_ID ', 'TIMESTEP_start ','AGE_start ','M1_start ','M2_start ','P_start ','LOGA_start ',&
                   'TIMESTEP_end ','AGE_end ','M1_end ','M2_end ','P_end ','LOGA_end ',&
                   'avgDM1R_msolpyr ','avgDM2R_msolpyr ',&
                   'avgDM1W_msolpyr ', 'avgDM2W_msolpyr ','CASE '
    121     format(2i8, 5E16.7, i8, 9E16.7, a20)

    open(114, file="./eve_txt/"//trim(zname)//"_DEATH.dat",&
              status='unknown', &
              access='sequential',&
              action='write')

    115     format(i8, 59E16.7)

    do while(err_status.eq.0)
        MODEL_ID = MODEL_ID + 1
        read(200, '(A)',iostat=err_status) fname_sneplot
        if(err_status.eq.-1) then
          ! if end of file, don't increment model_id and quit loop
          MODEL_ID = MODEL_ID -1
          exit
        endif

        write(100, 101) MODEL_ID, fname_sneplot
        !if(DEBUG_MODE.eq.2) write(*,*) "---READING SNEPLOT--- ", fname_sneplot

        read(200,*,iostat=err_status) modelimf, modeltype  ! read in IMF value and modeltype
        ! modeltype:
        !     -1 - single star
        !      0 - merger,
        !      1 - normal primary,
        !      2 - normal secondary,
        !      3 - single model,
        !      4 - binary hmg to get black hole mass too

        ! dealing with the weird format of the input files
        if(modeltype.ge.2) then ! Normal secondary file - has an extra line in the input file
            read(200,*,iostat=err_status) mixedimf, mixedage ! extra secondary information
            modelimf=modelimf-mixedimf ! calculate the IMF of non-rejuvenated and rejuvenated stars
            ! DEBUG ! print *,modelimf,mixedage
        endif
        if(modeltype.eq.4) then ! catch the QHE stars - extra initial binary orbit information
            read(200,*,iostat=err_status) initblackhole, initperiod ! so can put the black hole in at the end
            ! DEBUG ! print *, initblackhole, initperiod
        endif

        if(modelimf.lt.0d0) then
            modelimf=0d0 ! make sure we don't get -ve stars in IMF with rounding errors
        endif

        ! does the file we want to load exist? Put result in filecheck variable
        inquire(FILE=trim(models_dirname)//fname_sneplot,EXIST=filecheck)

        if(filecheck)then
            ! open sneplot (primary)
            open(210, FILE=trim(models_dirname)//fname_sneplot,&
                      STATUS='old',&
                      ACCESS='SEQUENTIAL',&
                      ACTION='read')

            ! Load the first line of the file to set up some variables
            read(210, 211, iostat=err_status_sneplot) (dummy(j),j=1,59)

            !if(DEBUG_MODE.eq.2) write(*,*) trim(fname_sneplot), modeltype, dummy(1), dummy(6)
            AGE_previous = dummy(2)
            m1_ZAMS = dummy(6)
            m2_ZAMS = dummy(38)  ! HFS: ZAMS mass of the secondary
            period_ZAMS = dummy(35) ! stays in years

            MT_bool = 0
            CEE_bool = 0
            MT_bool_summary = 0
            CEE_bool_summary = 0
            He_core_bool = .false.
            CO_core_bool = .false.

            do ! LOOP OVER ALL THE LINES IN THE MODEL TO LOCATE KEY EVOLUTION STEPS
                read(210, 211, iostat=err_status_sneplot) (dummy(j),j=1,59)

                ! MAIN SEQUENCE LIFETIME AND CASES
                if (dummy(7).gt.0.1) He_core_bool = .true.
                if (dummy(8).gt.0.1) CO_core_bool = .true.

                if (.not.He_core_bool) then
                  lifetime_MS = dummy(2)
                  case_str = 'A'
                endif
                if (He_core_bool.and..not.CO_core_bool) case_str = 'B'
                !if (He_core_bool.and.(CO_core_bool) case_str = 'C'
                if (CO_core_bool) case_str = 'C' ! HFS - @MAX will this work with your edge cases ?

                !write(*,*) case_str
                !if(DEBUG_MODE.eq.2) write(*,*) "CASE: ", case_str

                ! BINARY INTERACTIONS ONLY USEFUL FOR BINARIES!!!
                if(modeltype.ge.0.or.modeltype.le.2.or.modeltype.eq.4) then
                  ! QHE can interact (i think)
                     !!!!! STABLE MASS TRANSFER  !!!!
                     if((ABS(dummy(44)).gt.(1d-15)).and.(dummy(3).lt.dummy(36))) then
                         MT_bool_summary=1
                         if(CEE_bool.eq.1)then ! this will happend if stbale mass transfer after a CEE episode
                             ! record the CEE phase before moving on to the MT phase
                              CEE_bool=0 ! tell the code we out of CEE now
                              AGE_end = (dummy(2)+AGE_previous)/2 !the correct edge is between the two bins
                              avgDM1R_msolpyr = avgDM1R_msolpyr/rlfo_rows !(AGE_end-AGE_start)  
                              avgDM2R_msolpyr = avgDM2R_msolpyr/rlfo_rows !*(AGE_end-AGE_start) 
                              avgDM1W_msolpyr = avgDM1W_msolpyr/rlfo_rows !*(AGE_end-AGE_start) 
                              avgDM2W_msolpyr = avgDM2W_msolpyr/rlfo_rows !*(AGE_end-AGE_start) 


                              write(120, 121) MODEL_ID,&
                                               TIMESTEP_start, AGE_start, M1_start, M2_start, P_start, LOGA_start,&
                                               TIMESTEP_end, AGE_end, M1_end, M2_end, P_end, LOGA_end,&
                                               avgDM1R_msolpyr, avgDM2R_msolpyr,  avgDM1W_msolpyr, avgDM2W_msolpyr,&
                                               case_str


                         endif

                         if(MT_bool.eq.0) then ! first step of stable mass transfer
                             rlfo_rows = 1
                             TIMESTEP_start = dummy(1)
                             AGE_start = (dummy(2)+AGE_previous)/2
                             M1_start = dummy(6)
                             M2_start = dummy(38)
                             P_start = dummy(35)
                             LOGA_start = dummy(36)
                             ! NOT THE 1.989 FACTOR THAT WE THOUGHT WE NEEDED TO ADD!
                             avgDM1R_msolpyr = dummy(44)*3600*24*365.25
                             avgDM2R_msolpyr = dummy(45)*3600*24*365.25
                             avgDM1W_msolpyr = dummy(40)*3600*24*365.25
                             avgDM2W_msolpyr = dummy(41)*3600*24*365.25
                             MT_bool = 1
                             start_rlof_case=case_str

                         else if(MT_bool.eq.1) then ! continuing stable mass transfer
                             rlfo_rows = rlfo_rows + 1
                             TIMESTEP_end = dummy(1)
                             AGE_end = dummy(2)
                             M1_end = dummy(6)
                             M2_end = dummy(38)
                             P_end = dummy(35)
                             LOGA_end = dummy(36)
                             avgDM1R_msolpyr =avgDM1R_msolpyr+ dummy(44)*3600*24*365.25
                             avgDM2R_msolpyr =avgDM2R_msolpyr+ dummy(45)*3600*24*365.25
                             avgDM1W_msolpyr =avgDM1W_msolpyr+ dummy(40)*3600*24*365.25
                             avgDM2W_msolpyr =avgDM2W_msolpyr+ dummy(41)*3600*24*365.25
                        endif

                     else if(MT_bool.eq.1) then ! Stable mass transfer phase is over but DID happen
                         MT_bool=0 ! tell the code stable mass transfer is over
                         AGE_end = (dummy(2)+AGE_previous)/2 !the correct edge is between the two bins
                         avgDM1R_msolpyr = avgDM1R_msolpyr/rlfo_rows !(AGE_end-AGE_start)  
                         avgDM2R_msolpyr = avgDM2R_msolpyr/rlfo_rows !*(AGE_end-AGE_start) 
                         avgDM1W_msolpyr = avgDM1W_msolpyr/rlfo_rows !*(AGE_end-AGE_start) 
                         avgDM2W_msolpyr = avgDM2W_msolpyr/rlfo_rows !*(AGE_end-AGE_start) 

                         write(122, 121) MODEL_ID,&
                                         TIMESTEP_start, AGE_start, M1_start, M2_start, P_start, LOGA_start,&
                                         TIMESTEP_end, AGE_end, M1_end, M2_end, P_end, LOGA_end,&
                                         avgDM1R_msolpyr, avgDM2R_msolpyr,&
                                         avgDM1W_msolpyr, avgDM2W_msolpyr,start_rlof_case
                     endif ! end if Stable Mass transfer


                     !!!!    COMMON ENVELOPE    !!!!!

                     if(dummy(3).gt.dummy(36)) then
                         CEE_bool_summary=1
                         if(MT_bool.eq.1.and.DEBUG_MODE.eq.2) write(*,*) "PROBLEM - MT FLAG ON WHEN CEE"

                         if(CEE_bool.eq.0) then
                              rlfo_rows = 1
                              TIMESTEP_start = dummy(1)
                              AGE_start = (dummy(2)+AGE_previous)/2 ! the correct edge is between the 2 bins
                              M1_start = dummy(6)
                              M2_start = dummy(38)
                              P_start = dummy(35)
                              LOGA_start = dummy(36)
                              avgDM1R_msolpyr = dummy(44)*3600*24*365.25
                              avgDM2R_msolpyr = dummy(45)*3600*24*365.25
                              avgDM1W_msolpyr = dummy(40)*3600*24*365.25
                              avgDM2W_msolpyr = dummy(41)*3600*24*365.25
                              CEE_bool = 1
                              start_rlof_case = case_str
                         else if(CEE_bool.eq.1) then ! continuing stable mass transfer
                              rlfo_rows = rlfo_rows + 1
                              TIMESTEP_end = dummy(1)
                              AGE_end = dummy(2)
                              M1_end = dummy(6)
                              M2_end = dummy(38)
                              P_end = dummy(35)
                              LOGA_end = dummy(36)
                              avgDM1R_msolpyr = avgDM1R_msolpyr + dummy(44)*3600*24*365.25
                              avgDM2R_msolpyr = avgDM2R_msolpyr + dummy(45)*3600*24*365.25
                              avgDM1W_msolpyr = avgDM1W_msolpyr + dummy(40)*3600*24*365.25
                              avgDM2W_msolpyr = avgDM2W_msolpyr + dummy(41)*3600*24*365.25


                          endif

                     else if(CEE_bool.eq.1) then ! Stable mass transfer phase is over but DID happen
                          CEE_bool=0
                          AGE_end = (dummy(2)+AGE_previous)/2 !the correct edge is between the two bins
                          avgDM1R_msolpyr = avgDM1R_msolpyr/rlfo_rows!*(AGE_end-AGE_start)  
                          avgDM2R_msolpyr = avgDM2R_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 
                          avgDM1W_msolpyr = avgDM1W_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 
                          avgDM2W_msolpyr = avgDM2W_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 


                          write(120, 121) MODEL_ID,&
                                          TIMESTEP_start, AGE_start, M1_start, M2_start, P_start, LOGA_start,&
                                          TIMESTEP_end, AGE_end, M1_end, M2_end, P_end, LOGA_end,&
                                          avgDM1R_msolpyr, avgDM2R_msolpyr,&
                                          avgDM1W_msolpyr, avgDM2W_msolpyr,start_rlof_case

                     endif !! end ceeif
                     AGE_previous = dummy(2) 
                 endif !! end checking type if binary model that needs to look at interactions

                 if(err_status_sneplot.eq.-1.or.dummy(2).gt.1d11)then
                     ! end of file
                     ! recording the last row
                     write(114, 115) MODEL_ID, dummy
                     ! resolving any MT or CEE phases that have not been recorded
                     if (MT_bool.eq.1) then
                         MT_bool=0
                         AGE_end = (dummy(2)+AGE_previous)/2 !the correct edge is between the two bins
                         avgDM1R_msolpyr = avgDM1R_msolpyr/rlfo_rows!*(AGE_end-AGE_start)  
                         avgDM2R_msolpyr = avgDM2R_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 
                         avgDM1W_msolpyr = avgDM1W_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 
                         avgDM2W_msolpyr = avgDM2W_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 

                         write(122, 121) MODEL_ID,&
                                         TIMESTEP_start, AGE_start, M1_start, M2_start, P_start, LOGA_start,&
                                         TIMESTEP_end, AGE_end, M1_end, M2_end, P_end, LOGA_end,&
                                         avgDM1R_msolpyr, avgDM2R_msolpyr,&
                                         avgDM1W_msolpyr, avgDM2W_msolpyr,start_rlof_case

                     endif

                     if(CEE_bool.eq.1)then
                         CEE_bool=0
                         AGE_end = (dummy(2)+AGE_previous)/2 !the correct edge is between the two bins
                         avgDM1R_msolpyr = avgDM1R_msolpyr/rlfo_rows!*(AGE_end-AGE_start)  
                         avgDM2R_msolpyr = avgDM2R_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 
                         avgDM1W_msolpyr = avgDM1W_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 
                         avgDM2W_msolpyr = avgDM2W_msolpyr/rlfo_rows!*(AGE_end-AGE_start) 


                         write(120, 121) MODEL_ID,&
                                         TIMESTEP_start, AGE_start, M1_start, M2_start, P_start, LOGA_start,&
                                         TIMESTEP_end, AGE_end, M1_end, M2_end, P_end, LOGA_end,&
                                         avgDM1R_msolpyr, avgDM2R_msolpyr,&
                                         avgDM1W_msolpyr, avgDM2W_msolpyr,start_rlof_case
                                     
                    endif
                    ! Write SUMMARY table before we exit the loop

                    write(110, 111) MODEL_ID,&
                                  m1_ZAMS,m2_ZAMS,period_ZAMS,lifetime_MS,MT_bool_summary,CEE_bool_summary,&
                                  dummy(2), dummy(6), dummy(38), dummy(35),modelimf,mixedimf,modeltype

                    exit  ! if EOF exit
                 endif
            end do !File is over end, looping ober all the lines in sneplot
            ! Write to DEATH table and SUMMARY table before we go back to top of the loop
       endif ! end reading the sneplot file
       close(210) ! close the sneplot file

       ! IF CEE AT DEATH FLAG!

       ! THE LAST LINE CAN BE OUTPUT TO THE "DEATH" TABLE

       ! check SN type



    end do ! finish reading the input file and going through the models
end do!end of looking through tthe different metallicties

end program eve
