From 4d39b624b645399cf87969ebe1e78dbf89ac37f1 Mon Sep 17 00:00:00 2001 From: Taygun Bulmus Date: Sun, 23 Oct 2022 15:50:11 +0300 Subject: [PATCH] v1.1 --- .gitignore | 9 +- .vscode/settings.json | 2 + CHANGELOG | 13 + LICENSE | 2 +- collNuPy.py | 15 +- modules/RHS.py | 110 +++++++- modules/diffSolvers.py | 234 +++++++++++------- modules/hamilt.py | 8 +- modules/initial.py | 10 +- modules/plotgraphs.py | 28 ++- modules/readFiles.py | 5 +- modules/writeFiles.py | 18 +- parameters/physicalParameters.dat | 26 +- parameters/technicalParameters.dat | 11 +- scripts/bashScripts/plotgraphsScript.sh | 54 ---- scripts/create_SN_tXsdat.py | 130 ++++++++++ scripts/plotGraphsAfterSimulation.py | 79 ++++++ .../run_random_ri_magPolDecDist.py | 137 ++++++++++ scripts/run_random_ri_magPolDecDist.py | 103 ++++++++ 19 files changed, 798 insertions(+), 196 deletions(-) create mode 100644 CHANGELOG delete mode 100755 scripts/bashScripts/plotgraphsScript.sh create mode 100644 scripts/create_SN_tXsdat.py create mode 100644 scripts/plotGraphsAfterSimulation.py create mode 100644 scripts/random_ri_magPolDecDist/run_random_ri_magPolDecDist.py create mode 100644 scripts/run_random_ri_magPolDecDist.py diff --git a/.gitignore b/.gitignore index 47ed8d2..b588841 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,9 @@ +# Python Caches +*.pyc +*__pycache__*/ + +# Folders to ignore +backgroundProfiles/ +tmp/ results/ -parameters/old \ No newline at end of file +snSimulations/ \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json index 71b882e..b2c1584 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -4,12 +4,14 @@ ], "cSpell.ignoreWords": [ "Abbar", + "Athar", "DISTDIAG", "Emod", "Keil", "LSODA", "Termb", "Valeig", + "Xsdat", "cdist", "ebar", "figsize", diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 index 0000000..fd06a27 --- /dev/null +++ b/CHANGELOG @@ -0,0 +1,13 @@ +# v1.1 +- Updated gitignore. +- Converted `plotgraphsScript.sh` to python script (`plotGraphsAfterSimulation.py`). Remove `plotgraphsScript.sh`. +- Added scripts; + 1. SN post baryon density generator script from preSN file with added shockwave (`create_SN_tXsdat.py`). + 2. Run script for random `ri_km` and `magneticField_polynomialDecayDistance` (`run_random_ri_magPolDecDist.py`). +- Added `holdIntermediateData` option to technical parameters `technicalParameters.dat`. + - Added new functions without holding intermediate data `RHS.py`. + - Adopted `diffSolvers.py` for new `holdIntermediateData` variable. +- Fixed in matter hamiltonian (wrong assignment of Ye and nb when reading data from a file.) +- Minor changes in `plotgraphs.py`. +- Add function to copy background profile to results folder (`writeFiles.py`, `collnuPy.py`). +- Change all "STATUS" strings/variables to "INFO" in all code. \ No newline at end of file diff --git a/LICENSE b/LICENSE index f288702..a5bf218 100644 --- a/LICENSE +++ b/LICENSE @@ -300,7 +300,7 @@ or household purposes, or (2) anything designed or sold for incorporation into a dwelling. In determining whether a product is a consumer product, doubtful cases shall be resolved in favor of coverage. For a particular product received by a particular user, "normally used" refers to a -typical or common use of that class of product, regardless of the status +typical or common use of that class of product, regardless of the INFO of the particular user or of the way in which the particular user actually uses, or expects or is expected to use, the product. A product is a consumer product regardless of whether the product has substantial diff --git a/collNuPy.py b/collNuPy.py index b28e531..6d22d4a 100644 --- a/collNuPy.py +++ b/collNuPy.py @@ -33,7 +33,7 @@ # ============================ # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -51,7 +51,7 @@ class tLog: osCommand.makedirs(RESULTS_DIR) print(tLog.OK+"Results Folder Created.") else: - print(tLog.STATUS+"Results Folder Already Exists.") + print(tLog.INFO+"Results Folder Already Exists.") # ============================ # ============================ @@ -95,14 +95,14 @@ class tLog: # ============================ # Save final human readable data - # TODO: Comment out print after wrote the functions + # TODO: Comment out print after you wrote the functions if init.technicalParametersDic['output_humanReadable']: DATA_FOLDER_TXT= RESULTS_SIMULATION_DIR+ 'data/txt/' if not osCommand.path.exists(DATA_FOLDER_TXT): osCommand.makedirs(DATA_FOLDER_TXT) print(tLog.OK+"TXT Data Folder Created.") else: - print(tLog.STATUS+"TXT Data Folder Already Exists.") + print(tLog.INFO+"TXT Data Folder Already Exists.") # ============================ # rhoFlavAll.txt write2file.rhoFlavAll_txt(init,DATA_FOLDER_TXT+'rhoFlavAll.txt', rhoFlavAll) @@ -121,6 +121,11 @@ class tLog: #print(tLog.OK+"Saved: eigVal_eigVec.txt.") # ============================ + # ============================ + # Copy background profile file(s) + write2file.saveBackgroundProfiles(init, COLLECTIVE_NU_OSC_DIR+'backgroundProfiles/', DATA_FOLDER_NPZ) + # ============================ + # ============================ # Save final plots if init.technicalParametersDic['plotGraphs']: @@ -142,6 +147,6 @@ class tLog: osCommand.system('find . | grep -E "(__pycache__|.pyc|.pyo$)" | xargs rm -rf') print(tLog.OK+"Unused Files and Folders Removed.") #! ============================ - print(tLog.STATUS+'[FINISHED] =====> simulation'+str(simNum)+' <===== [FINISHED]\n') + print(tLog.INFO+'[FINISHED] =====> simulation'+str(simNum)+' <===== [FINISHED]\n') exit() # ============================ \ No newline at end of file diff --git a/modules/RHS.py b/modules/RHS.py index 6788845..cd56c7a 100644 --- a/modules/RHS.py +++ b/modules/RHS.py @@ -12,7 +12,7 @@ # ============================ # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -22,6 +22,20 @@ class tLog: # ============================ # RHS For Neutrinos f(r,rho) !!Only LSODA!! # Two Flavor +def rhs2Flav_noInterData(rhoAll, dist): + # ============================ + # Define the RHS + totHamAP = hamil.hamiltonian(rhoAll, dist).totalHam2Flav_km_1() + # rhoAll must be in the form of [2, energyMod, totFlav*2, totFlav*2] + + # Total Hamiltionian (totHam) [2, energyMod, totFlav, totFlav] + # For neutrinos -> totHam[0] + # For antineutrinos -> totHam[1] + + # Note: 3d array + 2d array gives a result that add all 2d array with each element of 3d array + # Take Commutator of density and Hamiltonian for each energy + return np.array([(-1j)* ((totHamAP[0] @ rhoAll[0])- (rhoAll[0] @ totHamAP[0]))\ + , (-1j)*((totHamAP[1] @ rhoAll[1])- (rhoAll[1] @ totHamAP[1]))]).flatten() def rhs2Flav(rhoAll, dist): # ============================ # Show intermediate distance and hold the values @@ -38,10 +52,10 @@ def rhs2Flav(rhoAll, dist): # To show current time currentTimeDate = datetime.datetime.now() fileObject = open('distanceAndTime.log','a+') - # Do not color the [STATUS] + # Do not color the [INFO] print('Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S"), file=fileObject) fileObject.close() - print(tLog.STATUS+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) + print(tLog.INFO+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) # ============================ # ============================ @@ -58,6 +72,20 @@ def rhs2Flav(rhoAll, dist): return np.array([(-1j)* ((totHamAP[0] @ rhoAll[0])- (rhoAll[0] @ totHamAP[0]))\ , (-1j)*((totHamAP[1] @ rhoAll[1])- (rhoAll[1] @ totHamAP[1]))]).flatten() # Three Flavor +def rhs3Flav_noInterData(rhoAll, dist): + # ============================ + # Define the RHS + totHamAP = hamil.hamiltonian(rhoAll, dist).totalHam3Flav_km_1() + # rhoAll must be in the form of [2, energyMod, totFlav*2, totFlav*2] + + # Total Hamiltionian (totHam) [2, energyMod, totFlav, totFlav] + # For neutrinos -> totHam[0] + # For antineutrinos -> totHam[1] + + # Note: 3d array + 2d array gives a result that add all 2d array with each element of 3d array + # Take Commutator of density and Hamiltonian for each energy + return np.array([(-1j)* ((totHamAP[0] @ rhoAll[0])- (rhoAll[0] @ totHamAP[0]))\ + , (-1j)*((totHamAP[1] @ rhoAll[1])- (rhoAll[1] @ totHamAP[1]))]).flatten() def rhs3Flav(rhoAll, dist): # ============================ # Show intermediate distance and hold the values @@ -74,10 +102,10 @@ def rhs3Flav(rhoAll, dist): # To show current time currentTimeDate = datetime.datetime.now() fileObject = open('distanceAndTime.log','a+') - # Do not color the [STATUS] + # Do not color the [INFO] print('Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S"), file=fileObject) fileObject.close() - print(tLog.STATUS+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) + print(tLog.INFO+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) # ============================ # ============================ @@ -94,6 +122,20 @@ def rhs3Flav(rhoAll, dist): return np.array([(-1j)* ((totHamAP[0] @ rhoAll[0])- (rhoAll[0] @ totHamAP[0]))\ , (-1j)*((totHamAP[1] @ rhoAll[1])- (rhoAll[1] @ totHamAP[1]))]).flatten() # Four Flavor +def rhs4Flav_noInterData(rhoAll, dist): + # ============================ + # Define the RHS + totHamAP = hamil.hamiltonian(rhoAll, dist).totalHam4Flav_km_1() + # rhoAll must be in the form of [2, energyMod, totFlav*2, totFlav*2] + + # Total Hamiltionian (totHam) [2, energyMod, totFlav, totFlav] + # For neutrinos -> totHam[0] + # For antineutrinos -> totHam[1] + + # Note: 3d array + 2d array gives a result that add all 2d array with each element of 3d array + # Take Commutator of density and Hamiltonian for each energy + return np.array([(-1j)* ((totHamAP[0] @ rhoAll[0])- (rhoAll[0] @ totHamAP[0]))\ + , (-1j)*((totHamAP[1] @ rhoAll[1])- (rhoAll[1] @ totHamAP[1]))]).flatten() def rhs4Flav(rhoAll, dist): # ============================ # Show intermediate distance and hold the values @@ -110,10 +152,10 @@ def rhs4Flav(rhoAll, dist): # To show current time currentTimeDate = datetime.datetime.now() fileObject = open('distanceAndTime.log','a+') - # Do not color the [STATUS] + # Do not color the [INFO] print('Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S"), file=fileObject) fileObject.close() - print(tLog.STATUS+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) + print(tLog.INFO+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) # ============================ # ============================ @@ -134,6 +176,20 @@ def rhs4Flav(rhoAll, dist): # ============================ # RHS For Neutrinos f(r,rho) !!Only LSODA!! # Two Flavor +def rhs2Flav_bigRho_noInterData(rhoAll, dist): + # ============================ + # Define the RHS + totHamAP = hamil.hamiltonian_BigRho(rhoAll, dist).totalHam2Flav_km_1() + # rhoAll must be in the form of [energyMod, totFlav*2, totFlav*2] + # Total Hamiltionian (totHam) [energyMod,totFlav*2,totFlav*2] + # totHam[N,0:totFlav,0:totFlav] => Neutrinos + # totHam[N,totFlav:totFlav*2,totFlav:totFlav*2] => Anti-Neutrinos + # totHam[N,0:totFlav,totFlav:totFlav*2] => Neutrino to Anti-neutrino + + # Note: 3d array + 2d array gives a result that add all 2d array with each element of 3d array + # Take Commutator of density and Hamiltonian for each energy + + return ((-1j)* ((totHamAP @ rhoAll)- (rhoAll @ totHamAP))).flatten() def rhs2Flav_bigRho(rhoAll, dist): # ============================ # Show intermediate distance and hold the values @@ -150,10 +206,10 @@ def rhs2Flav_bigRho(rhoAll, dist): # To show current time currentTimeDate = datetime.datetime.now() fileObject = open('distanceAndTime.log','a+') - # Do not color the [STATUS] + # Do not color the [INFO] print('Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S"), file=fileObject) fileObject.close() - print(tLog.STATUS+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) + print(tLog.INFO+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) # ============================ # ============================ @@ -170,6 +226,20 @@ def rhs2Flav_bigRho(rhoAll, dist): return ((-1j)* ((totHamAP @ rhoAll)- (rhoAll @ totHamAP))).flatten() # Three Flavor +def rhs3Flav_bigRho_noInterData(rhoAll, dist): + # ============================ + # Define the RHS + totHamAP = hamil.hamiltonian_BigRho(rhoAll, dist).totalHam3Flav_km_1() + # rhoAll must be in the form of [energyMod, totFlav*2, totFlav*2] + # Total Hamiltionian (totHam) [energyMod,totFlav*2,totFlav*2] + # totHam[N,0:totFlav,0:totFlav] => Neutrinos + # totHam[N,totFlav:totFlav*2,totFlav:totFlav*2] => Anti-Neutrinos + # totHam[N,0:totFlav,totFlav:totFlav*2] => Neutrino to Anti-neutrino + + # Note: 3d array + 2d array gives a result that add all 2d array with each element of 3d array + # Take Commutator of density and Hamiltonian for each energy + + return ((-1j)* ((totHamAP @ rhoAll)- (rhoAll @ totHamAP))).flatten() def rhs3Flav_bigRho(rhoAll, dist): # ============================ # Show intermediate distance and hold the values @@ -186,10 +256,10 @@ def rhs3Flav_bigRho(rhoAll, dist): # To show current time currentTimeDate = datetime.datetime.now() fileObject = open('distanceAndTime.log','a+') - # Do not color the [STATUS] + # Do not color the [INFO] print('Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S"), file=fileObject) fileObject.close() - print(tLog.STATUS+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) + print(tLog.INFO+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) # ============================ # ============================ @@ -206,6 +276,20 @@ def rhs3Flav_bigRho(rhoAll, dist): return ((-1j)* ((totHamAP @ rhoAll)- (rhoAll @ totHamAP))).flatten() # Four Flavor +def rhs4Flav_bigRho_noInterData(rhoAll, dist): + # ============================ + # Define the RHS + totHamAP = hamil.hamiltonian_BigRho(rhoAll, dist).totalHam4Flav_km_1() + # rhoAll must be in the form of [energyMod, totFlav*2, totFlav*2] + # Total Hamiltionian (totHam) [energyMod,totFlav*2,totFlav*2] + # totHam[N,0:totFlav,0:totFlav] => Neutrinos + # totHam[N,totFlav:totFlav*2,totFlav:totFlav*2] => Anti-Neutrinos + # totHam[N,0:totFlav,totFlav:totFlav*2] => Neutrino to Anti-neutrino + + # Note: 3d array + 2d array gives a result that add all 2d array with each element of 3d array + # Take Commutator of density and Hamiltonian for each energy + + return ((-1j)* ((totHamAP @ rhoAll)- (rhoAll @ totHamAP))).flatten() def rhs4Flav_bigRho(rhoAll, dist): # ============================ # Show intermediate distance and hold the values @@ -222,10 +306,10 @@ def rhs4Flav_bigRho(rhoAll, dist): # To show current time currentTimeDate = datetime.datetime.now() fileObject = open('distanceAndTime.log','a+') - # Do not color the [STATUS] + # Do not color the [INFO] print('Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S"), file=fileObject) fileObject.close() - print(tLog.STATUS+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) + print(tLog.INFO+'Current distance is %.7f' % dist, 'km and time :', currentTimeDate.strftime("%Y-%m-%d %H:%M:%S")) # ============================ # ============================ diff --git a/modules/diffSolvers.py b/modules/diffSolvers.py index f4f68fd..13704af 100644 --- a/modules/diffSolvers.py +++ b/modules/diffSolvers.py @@ -8,7 +8,6 @@ import time import datetime import os as osCommand -from distutils.dir_util import copy_tree from shutil import rmtree import numpy as np import RHS as func @@ -21,7 +20,7 @@ # ============================ # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -55,14 +54,15 @@ def odeintwSolver(init, RESULTS_SIMULATION_DIR): else: print(tLog.ERROR+"Data Folder Already Exists. Check The Folder: "+ DATA_FOLDER) print(tLog.EXIT) - # Create intermediateData folder - INTERMEDIATE_DATA_FOLDER= RESULTS_SIMULATION_DIR+ 'data/intermediateData/' - if not osCommand.path.exists(INTERMEDIATE_DATA_FOLDER): - osCommand.makedirs(INTERMEDIATE_DATA_FOLDER) - print(tLog.OK+"Intermediate Data Folder Created.") - else: - print(tLog.ERROR+"Intermediate Data Folder Already Exists. Check The Folder: "+ INTERMEDIATE_DATA_FOLDER) - print(tLog.EXIT) + if init.technicalParametersDic['holdIntermediateData']: + # Create intermediateData folder + INTERMEDIATE_DATA_FOLDER= RESULTS_SIMULATION_DIR+ 'data/intermediateData/' + if not osCommand.path.exists(INTERMEDIATE_DATA_FOLDER): + osCommand.makedirs(INTERMEDIATE_DATA_FOLDER) + print(tLog.OK+"Intermediate Data Folder Created.") + else: + print(tLog.ERROR+"Intermediate Data Folder Already Exists. Check The Folder: "+ INTERMEDIATE_DATA_FOLDER) + print(tLog.EXIT) # Change directory to data # Code saves intermediate data osCommand.chdir(DATA_FOLDER) @@ -75,71 +75,138 @@ def odeintwSolver(init, RESULTS_SIMULATION_DIR): # ============================ # Solve E.o.M. # Starting time - print(tLog.STATUS+'Start solving E.o.M with odeintw. Current date and time :', timeInitialDate.strftime("%Y-%m-%d %H:%M:%S")) - if init.dim_rho_2totFlav_Bool: - # Default error tolerances - if (init.technicalParametersDic['tolerance_relativeError'] == 0 and init.technicalParametersDic['tolerance_absoluteError'] == 0): - if init.totFlav == 2: - # Two Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_bigRho, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) - elif init.totFlav == 3: - # Three Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_bigRho, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) - elif init.totFlav == 4: - # Four Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_bigRho, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) - # User defined error tolerances + print(tLog.INFO+'Start solving E.o.M with odeintw. Current date and time :', timeInitialDate.strftime("%Y-%m-%d %H:%M:%S")) + # Hold intermediate data + if init.technicalParametersDic['holdIntermediateData']: + if init.dim_rho_2totFlav_Bool: + # Default error tolerances + if (init.technicalParametersDic['tolerance_relativeError'] == 0 and init.technicalParametersDic['tolerance_absoluteError'] == 0): + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_bigRho, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_bigRho, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_bigRho, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + # User defined error tolerances + else: + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_bigRho, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_bigRho, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_bigRho, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) else: - if init.totFlav == 2: - # Two Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_bigRho, init.rhoInit_flav\ - , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ - , atol=init.technicalParametersDic['tolerance_absoluteError']\ - , mxstep=500000000) - elif init.totFlav == 3: - # Three Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_bigRho, init.rhoInit_flav\ - , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ - , atol=init.technicalParametersDic['tolerance_absoluteError']\ - , mxstep=500000000) - elif init.totFlav == 4: - # Four Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_bigRho, init.rhoInit_flav\ - , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ - , atol=init.technicalParametersDic['tolerance_absoluteError']\ - , mxstep=500000000) + # Default error tolerances + if (init.technicalParametersDic['tolerance_relativeError'] == 0 and init.technicalParametersDic['tolerance_absoluteError'] == 0): + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + # User defined error tolerances + else: + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) else: - # Default error tolerances - if (init.technicalParametersDic['tolerance_relativeError'] == 0 and init.technicalParametersDic['tolerance_absoluteError'] == 0): - if init.totFlav == 2: - # Two Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) - elif init.totFlav == 3: - # Three Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) - elif init.totFlav == 4: - # Four Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) - # User defined error tolerances + if init.dim_rho_2totFlav_Bool: + # Default error tolerances + if (init.technicalParametersDic['tolerance_relativeError'] == 0 and init.technicalParametersDic['tolerance_absoluteError'] == 0): + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_bigRho_noInterData, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_bigRho_noInterData, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_bigRho_noInterData, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + # User defined error tolerances + else: + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_bigRho_noInterData, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_bigRho_noInterData, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_bigRho_noInterData, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) else: - if init.totFlav == 2: - # Two Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav, init.rhoInit_flav\ - , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ - , atol=init.technicalParametersDic['tolerance_absoluteError']\ - , mxstep=500000000) - elif init.totFlav == 3: - # Three Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav, init.rhoInit_flav\ - , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ - , atol=init.technicalParametersDic['tolerance_absoluteError']\ - , mxstep=500000000) - elif init.totFlav == 4: - # Four Flavor - rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav, init.rhoInit_flav\ - , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ - , atol=init.technicalParametersDic['tolerance_absoluteError']\ - , mxstep=500000000) + # Default error tolerances + if (init.technicalParametersDic['tolerance_relativeError'] == 0 and init.technicalParametersDic['tolerance_absoluteError'] == 0): + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_noInterData, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_noInterData, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_noInterData, init.rhoInit_flav, init.distAll_km, full_output=True, mxstep=500000000) + # User defined error tolerances + else: + if init.totFlav == 2: + # Two Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs2Flav_noInterData, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 3: + # Three Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs3Flav_noInterData, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) + elif init.totFlav == 4: + # Four Flavor + rhoFinalAll,infow = integralSolver.odeintw(func.rhs4Flav_noInterData, init.rhoInit_flav\ + , init.distAll_km, full_output=True, rtol=init.technicalParametersDic['tolerance_relativeError']\ + , atol=init.technicalParametersDic['tolerance_absoluteError']\ + , mxstep=500000000) # ============================ # ============================ @@ -161,25 +228,25 @@ def odeintwSolver(init, RESULTS_SIMULATION_DIR): minutes = tmpTimeCalculation // 60 tmpTimeCalculation %= 60 seconds = tmpTimeCalculation - print(tLog.STATUS+"Total Process Time, %6.2fs, d:h:m:s : %d:%d:%d:%d"% (tmpTimeCalculationSec, days, hours, minutes, seconds)) + print(tLog.INFO+"Total Process Time, %6.2fs, d:h:m:s : %d:%d:%d:%d"% (tmpTimeCalculationSec, days, hours, minutes, seconds)) # ============================ # Save Total Process Time to file - # Do not color the [STATUS] + # Do not color the [INFO] fileObject = open('distanceAndTime.log','a+') print('Total Process Time, %6.2fs, d:h:m:s : %d:%d:%d:%d'% (tmpTimeCalculationSec, days, hours, minutes, seconds), file=fileObject) fileObject.close() # ============================ # Number Of Evolution - print(tLog.STATUS+'Number of evaluation: %d' % infow['nfe'][-1]) + print(tLog.INFO+'Number of evaluation: %d' % infow['nfe'][-1]) # ============================ # Check differential equation solver method changed or not methodIndicators = infow['mused'] if methodIndicators[0] == 1: - print(tLog.STATUS+'Differential equation solver started with adams (nonstiff) method.') + print(tLog.INFO+'Differential equation solver started with adams (nonstiff) method.') else: - print(tLog.STATUS+'Differential equation solver started with bdf (stiff) method.') - indicatorChangePrint = [tLog.STATUS+'Differential equation solver method is adams (nonstiff)'\ - , tLog.STATUS+'Differential equation solver method is bdf (stiff) '] + print(tLog.INFO+'Differential equation solver started with bdf (stiff) method.') + indicatorChangePrint = [tLog.INFO+'Differential equation solver method is adams (nonstiff)'\ + , tLog.INFO+'Differential equation solver method is bdf (stiff) '] # ============================ # If tmpIndicatorCheckControl == 1, the method changed. tmpIndicatorCheckControl = 0 @@ -204,9 +271,9 @@ def odeintwSolver(init, RESULTS_SIMULATION_DIR): # ============================ # Print if the method is changed or not if tmpIndicatorCheckControl == 0: - print(tLog.STATUS+'Differential equation solver method have not changed.') + print(tLog.INFO+'Differential equation solver method have not changed.') elif tmpIndicatorCheckControl == 1: - print(tLog.STATUS+'Differential equation solver method have changed.') + print(tLog.INFO+'Differential equation solver method have changed.') # ============================ ''' ALL INFOS @@ -229,8 +296,9 @@ def odeintwSolver(init, RESULTS_SIMULATION_DIR): 1: adams (nonstiff), 2: bdf (stiff) ''' # ============================ - # Remove intermediateData folder - rmtree(INTERMEDIATE_DATA_FOLDER) + if init.technicalParametersDic['holdIntermediateData']: + # Remove intermediateData folder + rmtree(INTERMEDIATE_DATA_FOLDER) # Return Back To RESULTS_SIMULATION_DIR Folder osCommand.chdir(RESULTS_SIMULATION_DIR) diff --git a/modules/hamilt.py b/modules/hamilt.py index 5e76127..dd0369e 100644 --- a/modules/hamilt.py +++ b/modules/hamilt.py @@ -16,7 +16,7 @@ # ============================ # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -125,8 +125,9 @@ class tLog: # ============================ # Interpolate n_baryon and Y_e w.r.t red distance - nBaryon_1_km3_LAMBDA = interp1d(distRead_km, tmp0, kind='linear'); del tmp0 - Y_e_LAMBDA = interp1d(distRead_km, tmp1, kind='linear'); del tmp1 + Y_e_LAMBDA = interp1d(distRead_km, tmp0, kind='linear'); del tmp0 + nBaryon_1_km3_LAMBDA = interp1d(distRead_km, tmp1, kind='linear'); del tmp1 + else: print(tLog.ERROR+'The file '+SN_MATTER_PROFILE_DATA_DIR+'does not exist.') print(tLog.ERROR+'Please check the file name and path.') @@ -184,6 +185,7 @@ class tLog: # ============================ # ============================ +# If dimension is totFlav x totFlav class hamiltonian: # ============================ # Variables diff --git a/modules/initial.py b/modules/initial.py index 16e51b5..caeae68 100644 --- a/modules/initial.py +++ b/modules/initial.py @@ -15,7 +15,7 @@ # ============================ # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -58,7 +58,7 @@ class tLog: #! ============================ if not sterileBool and nuDistParam == 5: # ERROR Sterile neutrinos are not considered, but neutrino_distributionParameter=5 - print(tLog.STATUS+'Sterile neutrinos are not considered, but neutrino_distributionParameter=5.') + print(tLog.INFO+'Sterile neutrinos are not considered, but neutrino_distributionParameter=5.') exit(tLog.EXIT) elif sterileBool and nuDistParam != 5: # WARNING @@ -255,12 +255,12 @@ class tLog: # ERROR Temperature Errors if 0 in tempArray_MeV and nuDistParam != 5: if nuDistParam == 1: - print(tLog.STATUS+'Fermi Dirac Distribution: (E^2)/(F_2(0) T^3 (exp(E/T) +1)) <== [Duan:2006an] eqn (12)') + print(tLog.INFO+'Fermi Dirac Distribution: (E^2)/(F_2(0) T^3 (exp(E/T) +1)) <== [Duan:2006an] eqn (12)') print(tLog.ERROR+'Temperature of one of the neutrinos is zero. Any of temperatures can not be zero ' 'for Fermi Dirac Distribution.') exit(tLog.EXIT) elif nuDistParam == 2: - print(tLog.STATUS+'Pinched Dirac Distribution: (((alpha+1)/)^alpha)*(E^alpha/Gamma(alpha+1))*exp(-(alpha+1)*E/) [Keil:2002in]') + print(tLog.INFO+'Pinched Dirac Distribution: (((alpha+1)/)^alpha)*(E^alpha/Gamma(alpha+1))*exp(-(alpha+1)*E/) [Keil:2002in]') print(tLog.ERROR+'Temperature of one of the neutrinos is zero. Any of temperatures can not be zero ' 'for Pinched Dirac Distribution.') exit(tLog.EXIT) @@ -271,7 +271,7 @@ class tLog: or physicalParametersDic['temperature_eb'] == 0\ or physicalParametersDic['temperature_mub'] == 0\ or physicalParametersDic['temperature_taub'] == 0: - print(tLog.STATUS+'Fermi Dirac Distribution: (E^2)/(F_2(0) T^3 (exp(E/T) +1)) <== [Duan:2006an] eqn (12)') + print(tLog.INFO+'Fermi Dirac Distribution: (E^2)/(F_2(0) T^3 (exp(E/T) +1)) <== [Duan:2006an] eqn (12)') print(tLog.ERROR+'Temperature of one of the active neutrinos is zero. Any of temperatures can not be zero ' 'for Fermi Dirac Distribution.') exit(tLog.EXIT) diff --git a/modules/plotgraphs.py b/modules/plotgraphs.py index 79b6666..858255a 100644 --- a/modules/plotgraphs.py +++ b/modules/plotgraphs.py @@ -7,7 +7,7 @@ # ============================ # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -61,7 +61,7 @@ def __init__(self, SIMULATION_PATH): osCommand.makedirs(self.FIGURE_PATH) print(tLog.OK+"figures Folder Created.") else: - print(tLog.STATUS+"figures Folder Already Exists.") + print(tLog.INFO+"figures Folder Already Exists.") # ============================ # ============================ @@ -120,7 +120,7 @@ def distDiag(self): osCommand.makedirs(FIG_DISTDIAG_PATH) print(tLog.OK+"figures/distDiag/ Folder Created.") else: - print(tLog.STATUS+"figures/distDiag/ Folder Already Exists.") + print(tLog.INFO+"figures/distDiag/ Folder Already Exists.") # Create Y Axises yAxises = np.zeros((len(self.distAll_km), self.Emod, self.totFlav*2)) @@ -157,7 +157,10 @@ def distDiag(self): # Clear Figure plt.clf() - plt.close('all') + plt.close('all') + # ============================ + print(tLog.OK+"Diagonal elements of density matrix to distance graph(s) for each energy are created.") + # ============================ def energyDiag(self): # Paths FIG_DISTDIAG_PATH = self.FIGURE_PATH+'energyDiag/' @@ -165,7 +168,7 @@ def energyDiag(self): osCommand.makedirs(FIG_DISTDIAG_PATH) print(tLog.OK+"figures/energyDiag/ Folder Created.") else: - print(tLog.STATUS+"figures/energyDiag/ Folder Already Exists.") + print(tLog.INFO+"figures/energyDiag/ Folder Already Exists.") # Create Y Axises yAxisesBeg = np.zeros((self.Emod, self.totFlav*2)) @@ -214,6 +217,7 @@ def energyDiag(self): # Clear Figure plt.clf() plt.close('all') + print(tLog.OK+"Diagonal elements of final density matrix to energy graph(s) are created.") def distHamiltDiag(self): # ============================ # hamiltonians.npz load @@ -232,7 +236,7 @@ def distHamiltDiag(self): osCommand.makedirs(FIG_DISTDIAG_PATH) print(tLog.OK+"figures/distHamiltDiag/ Folder Created.") else: - print(tLog.STATUS+"figures/distHamiltDiag/ Folder Already Exists.") + print(tLog.INFO+"figures/distHamiltDiag/ Folder Already Exists.") # ============================ # ============================ @@ -251,13 +255,14 @@ def distHamiltDiag(self): , color='r'\ , alpha=0.4) # Labels - plt.ylabel('$H_{osc}(0,0)$ [km$^{-1}$]') + plt.ylabel('$H_{osc}(1,1)$ [km$^{-1}$]') plt.xlabel('Distance [km]') # Save The Graph plt.savefig(FIG_DISTDIAG_PATH+'distHamiltDiag_Osc.'+self.pltFormat) # Clear Figure plt.clf() plt.close('all') + print(tLog.OK+"(0,0) elements of oscillation Hamiltonian (km^-1) to distance graph(s) for each energy are created.") # ============================ # Matter Hamiltonian if hamBoolAll_OscMatEMSelfSA[1]: @@ -267,7 +272,7 @@ def distHamiltDiag(self): plt.plot(self.distAll_km, (np.abs(np.real(hamMat_km_1[:, 1, 1]))), color='r', label='$H_{mat}(2,2)$') if self.totFlav == 3: plt.plot(self.distAll_km, (np.abs(np.real(hamMat_km_1[:, 2, 2]))), color='b', label='$H_{mat}(3,3)$') - plt.ylabel('$H_{mat}(0,0)$ [km$^{-1}$]') + plt.ylabel('$H_{mat}$ [km$^{-1}$]') # Labels and legend plt.xlabel('Distance [km]') plt.legend(loc = 'upper right') @@ -276,6 +281,7 @@ def distHamiltDiag(self): # Clear Figure plt.clf() plt.close('all') + print(tLog.OK+"(0,0), (1,1) elements of matter Hamiltonian (km^-1) to distance graph(s) are created.") # ============================ # EM Hamiltonian if hamBoolAll_OscMatEMSelfSA[2]: @@ -291,13 +297,14 @@ def distHamiltDiag(self): # Clear Figure plt.clf() plt.close('all') + print(tLog.OK+"(0,-1) elements of electromagnetic Hamiltonian (km^-1) or muB to distance graph(s) are created.") # ============================ # SelfSA Hamiltonian if hamBoolAll_OscMatEMSelfSA[3]: hamSelfSA_km_1= hamiltonians_npz['hamSelfSA_km_1'] # Plot - plt.plot(self.distAll_km, (np.abs(np.real(hamSelfSA_km_1[:, 0, 0]))), color='k', label='$H_{mat}(0,0)$') - plt.ylabel('$H_{mat}(0,0)$ [km$^{-1}$]') + plt.plot(self.distAll_km, (np.abs(np.real(hamSelfSA_km_1[:, 0, 0]))), color='k', label='$H_{self}(1,1)$') + plt.ylabel('$H_{self}(1,1)$ [km$^{-1}$]') # Labels and legend plt.xlabel('Distance [km]') plt.legend(loc = 'upper right') @@ -306,4 +313,5 @@ def distHamiltDiag(self): # Clear Figure plt.clf() plt.close('all') + print(tLog.OK+"(0,0) elements of self interaction Hamiltonian (km^-1) to distance graph(s) are created.") # ============================ diff --git a/modules/readFiles.py b/modules/readFiles.py index 73979ea..44a1527 100644 --- a/modules/readFiles.py +++ b/modules/readFiles.py @@ -7,7 +7,7 @@ import numpy as np # Colors For Printing class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -36,6 +36,7 @@ def technicalParameters_read(): # With Number Comment # Same for two and three flavor technicalParametersNames= ('holdData_Every', + 'holdIntermediateData', 'plotGraphs', 'plotGraphs_diagonalRhoAllEnergy_2distance', 'plotGraphs_diagonalRhoFinal_2energy', @@ -49,7 +50,7 @@ def technicalParameters_read(): 'tolerance_absoluteError') # If parameters' value must be 0 or 1, add their indexes below - err01=('plotGraphs','plotGraphs_diagonalRhoAllEnergy_2distance','plotGraphs_diagonalRhoFinal_2energy'\ + err01=('holdIntermediateData','plotGraphs','plotGraphs_diagonalRhoAllEnergy_2distance','plotGraphs_diagonalRhoFinal_2energy'\ ,'plotGraphs_hamiltonianAllEnergy_2distance'\ ,'output_humanReadable','output_distanceHamiltonianAllE', 'output_distance_eigenValuesAllE') diff --git a/modules/writeFiles.py b/modules/writeFiles.py index ec87c69..de4b532 100644 --- a/modules/writeFiles.py +++ b/modules/writeFiles.py @@ -6,10 +6,11 @@ """ import numpy as np import hamilt as hamil +from shutil import copy2 # ============================ class tLog: - STATUS = '\033[94m'+'[STATUS] '+'\033[0m' + INFO = '\033[94m'+'[INFO] '+'\033[0m' OK = '\033[92m'+'[OK] '+'\033[0m' WARNING = '\033[93m'+'[WARNING] '+'\033[0m' ERROR = '\033[91m'+'[ERROR] '+'\033[0m' @@ -275,3 +276,18 @@ def eigVal_eigVec_txt(init, fullFilePath, rhoFlavAll): print(tLog.ERROR+'eigVal_eigVec.txt file not implemented yet.') # ============================ +# ============================ +# Save Background Profile +def saveBackgroundProfiles(init, backgroundFolderPath, saveFolderPath): + # Do you use default background profile? + if not init.physicalParametersDic['use_defaultMatterProfile']: + matterProfile_fileName= str(init.physicalParametersDic['matterProfile_fileName'][:-1]) + # Save file + copy2(backgroundFolderPath+matterProfile_fileName, saveFolderPath) + print(tLog.INFO+'Matter background profile file is copied to: '+saveFolderPath) + if not init.physicalParametersDic['use_defaultMagneticProfile']: + magneticProfile_fileName= str(init.physicalParametersDic['magneticProfile_fileName'][:-1]) + # Save file + copy2(backgroundFolderPath+magneticProfile_fileName, saveFolderPath) + print(tLog.INFO+'EM background profile file is copied to: '+saveFolderPath) +# ============================ \ No newline at end of file diff --git a/parameters/physicalParameters.dat b/parameters/physicalParameters.dat index 63e84b6..80dedfa 100644 --- a/parameters/physicalParameters.dat +++ b/parameters/physicalParameters.dat @@ -11,12 +11,12 @@ hamiltonian_electromagnetic=1 # Valid only (0 || 1) # ------------------------------------------------------------------------------ flavor_number=2 # Valid only (2, 3, 4) => 2(e,x), 3(e,mu,tau), 4(e,mu,tau,sterile) -neutrino_hierarchy=1 # Valid only (1 normal, -1 inverted) +neutrino_hierarchy=1 # Valid only (1 normal, -1 inverted) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ theta12=0.146 # [Radian] Valid for (flavor_number=2,3,4) , (hamiltonian_oscillation=1) -theta23=0.146 # [Radian] Valid for (flavor_number=3,4) , (hamiltonian_oscillation=1) +theta23=0 # [Radian] Valid for (flavor_number=3,4) , (hamiltonian_oscillation=1) theta13=0 # [Radian] Valid for (flavor_number=3,4) , (hamiltonian_oscillation=1) theta_14=0 # [Radian] Valid for (flavor_number=4) , (hamiltonian_oscillation=1) theta_24=0 # [Radian] Valid for (flavor_number=4) , (hamiltonian_oscillation=1) @@ -26,20 +26,20 @@ CP_Phase=0 # [Radian] , (hamiltonian_oscillation=1) # ------------------------------------------------------------------------------ deltaM_square21=2.56e-15 # [MeV^2] Valid for (flavor_number=2,3,4) -deltaM_square32=2.40e-13 # [MeV^2] Valid for (flavor_number=3,4) +deltaM_square32=0 # [MeV^2] Valid for (flavor_number=3,4) deltaM_square41=0 # [MeV^2] Valid for (flavor_number=4) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ -neutrino_MagneticMoment=1e-15 # [Bohr Magneton] Valid for (hamiltonian_electromagnetic=1) +neutrino_MagneticMoment=5e-15 # [Bohr Magneton] Valid for (hamiltonian_electromagnetic=1) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ -numberOf_energyMode=1 # Valid for only integer -energy_initial=1 # [MeV] -energy_final=1 # [MeV] -distance_initial=50.03995115698231 # [km] -distance_final=3200 # [km] +numberOf_energyMode=1 # Valid for only integer +energy_initial=1.0 # [MeV] +energy_final=1.0 # [MeV] +distance_initial=49.972151143341506 # [km] +distance_final=54 # [km] neutrino_distributionParameter=7 # Valid only (1,2,3,4,5,6,7), For profiles see below # ------------------------------------------------------------------------------ @@ -72,14 +72,14 @@ temperature_sterileb=0 # [MeV] Valid for (flavor_number=4, neutrino_densityParam # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ -use_defaultMatterProfile=1 # Valid for (0 || 1), (hamiltonian_matter=1) +use_defaultMatterProfile=0 # Valid for (0 || 1), (hamiltonian_matter=1) matterDensity_profile=1 # Valid for only (0,1,2) => (Constant, exponential, polynomial) (use_defaultMatterProfile=1) matterDensity_initial=1.8e+7 # [g/cm^3] Valid for (use_defaultMatterProfile=1) matterDensity_exponentialDecay=200 # [km] Valid for (use_defaultMatterProfile=1, matterDensity_profile=1) matterDensity_polynomialDecay=3 # Valid for (use_defaultMatterProfile=1, matterDensity_profile=2) electronFraction_constant=0.45 # Valid for (use_defaultMatterProfile=1) -matterProfile_fileName=exampleProfile.dat # Valid for (use_defaultMatterProfile=0) It must be under "backgroundProfiles" folder. +matterProfile_fileName=SN_distCM_Ye_Nb_g_cm3_t5.0s.dat # Valid for (use_defaultMatterProfile=0) It must be under "backgroundProfiles" folder. # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ @@ -88,7 +88,7 @@ magneticField_profile=2 # Valid for only (0,1,2) => (Constant, magneticField_initial=1e+15 # [Gauss] Valid for (hamiltonian_electromagnetic=1, use_defaultMagneticProfile=1) magneticField_exponentialDecay=200 # [km] Valid for (use_defaultMagneticProfile=1, magneticField_profile=1) magneticField_polynomialDecayPower=2 # Valid for (use_defaultMagneticProfile=1, magneticField_profile=2) -magneticField_polynomialDecayDistance=49.977992853062304 # [km] Valid for (use_defaultMagneticProfile=1, magneticField_profile=2) +magneticField_polynomialDecayDistance=49.9954430063595 # [km] Valid for (use_defaultMagneticProfile=1, magneticField_profile=2) magneticField_fileName=magneticField_profile.dat # Valid for (use_defaultMagneticProfile=0). # ------------------------------------------------------------------------------ @@ -124,4 +124,4 @@ couplingConstant_sterile=1 # [Fermi_Constant] Valid for (flavor_number=4, hamilt # [6] Only eBox [] : Electron neutrino box spectrum # [7] Only ebarBox [] : Anti Electron neutrino box spectrum # ----------------------------------------------------------------------------- -# ============================================================================== \ No newline at end of file +# ============================================================================== diff --git a/parameters/technicalParameters.dat b/parameters/technicalParameters.dat index 9019686..6c058ff 100644 --- a/parameters/technicalParameters.dat +++ b/parameters/technicalParameters.dat @@ -1,24 +1,25 @@ # ------------------------------------------------------------------------------ -holdData_Every=1 # [km] +holdData_Every=10 # [km] +holdIntermediateData=1 # Valid for (0 || 1) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ plotGraphs= 1 # Valid only (0 || 1) plotGraphs_diagonalRhoAllEnergy_2distance=1 # Valid for (0 || 1), (plotGraphs= 1) -plotGraphs_diagonalRhoFinal_2energy=1 # Valid for (0 || 1), (plotGraphs= 1) -plotGraphs_hamiltonianAllEnergy_2distance=1 # Valid for (0 || 1), (plotGraphs= 1) +plotGraphs_diagonalRhoFinal_2energy=0 # Valid for (0 || 1), (plotGraphs= 1) +plotGraphs_hamiltonianAllEnergy_2distance=0 # Valid for (0 || 1), (plotGraphs= 1) plot_savingFormat=png # Valid for (png,rgba,eps,svgz,pgf,ps,raw,svg,pdf), (plotGraphs= 1) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ output_distanceHamiltonianAllE=1 # Valid only (0 || 1) -output_distance_eigenValuesAllE=1 # Valid only (0 || 1) +output_distance_eigenValuesAllE=0 # Valid only (0 || 1) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ -output_humanReadable=1 # Valid only (0 || 1) +output_humanReadable=0 # Valid only (0 || 1) # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ diff --git a/scripts/bashScripts/plotgraphsScript.sh b/scripts/bashScripts/plotgraphsScript.sh deleted file mode 100755 index ac445bd..0000000 --- a/scripts/bashScripts/plotgraphsScript.sh +++ /dev/null @@ -1,54 +0,0 @@ -#!/usr/bin/env bash - -cd ../../ -# collnu Directory -CURR_DIR=$(pwd) -if [ -z "$1" ] - then - read -p "Enter simulation name: " SIMU_NAME - else - SIMU_NAME=$1 -fi - -RESULTS_DIR = "$CURR_DIR/results/" -NPZ_FILE_PATH="$CURR_DIR/results/$SIMU_NAME/data/rhoFlavAll.npz" - -echo "Try to reach file destination:" -echo $NPZ_FILE_PATH - -### Check if a directory does not exist ### -if [ ! -f "$NPZ_FILE_PATH" ]; then - echo "Simulation $SIMU_NAME DOES NOT exists." - exit 9999 # die with error code 9999 -fi -echo "Reached file destination:" - -echo "Avaible plots are following" -echo "Diagonal elements of density matrix to energy => plot_EnergyDiag : a" -echo "Diagonal elements of density matrix to distance plot => plot_DistDiag : b" -echo "" -echo "Example: To plot plot_EnergyDiag and plot_Adiabaticity, enter ag" - -read -p "Which graph do you want to plot? " plotCheck - -python3 <<- EOF -from sys import path -path.insert(0, "$CURR_DIR/modules") -import plotgraphs as plot -print("The plotting is started.") - -RESULTS_SIMULATION_DIR= "$CURR_DIR/results/$SIMU_NAME/" - -# ====================================================================== -# PLOTS -# ====================================================================== -if 'a' in "$plotCheck": - print('RUN energyDiag') - plot.plotMain(RESULTS_SIMULATION_DIR).energyDiag() - print('Finished plot_energyDiag') -if 'b' in "$plotCheck": - print('RUN distDiag') - plot.plotMain(RESULTS_SIMULATION_DIR).distDiag() - print('Finished plot_distDiag') -# ====================================================================== -EOF diff --git a/scripts/create_SN_tXsdat.py b/scripts/create_SN_tXsdat.py new file mode 100644 index 0000000..e5aa857 --- /dev/null +++ b/scripts/create_SN_tXsdat.py @@ -0,0 +1,130 @@ +import numpy as np +import matplotlib.pyplot as plt +import argparse +from pathlib import Path +import os as osCommand # File operations + +# ============================ +# Colors For Printing +class tLog: + INFO = '\033[94m'+'[INFO] '+'\033[0m' + OK = '\033[92m'+'[OK] '+'\033[0m' + WARNING = '\033[93m'+'[WARNING] '+'\033[0m' + ERROR = '\033[91m'+'[ERROR] '+'\033[0m' + EXIT = '\033[91m'+'[EXIT] create_SN_tXsdat.py' + '\033[0m' +# ============================ + +# ============================ +# Create Ye as a Step Function +def stepYe_for_Athar1995cx(distAll_km, innerStep_km=695.700, outerStep_km=695700\ + , innerYe=0.35, midYe=0.4998, outerYe=0.75): + print(tLog.INFO+'Creating Ye as a Step Function.') + Ye= np.zeros(len(distAll_km)) + for i1 in range(len(distAll_km)): + if (distAll_km[i1] <= innerStep_km): + Ye[i1]= innerYe + elif (distAll_km[i1] <= outerStep_km and distAll_km[i1] > innerStep_km): + Ye[i1]= midYe + elif (distAll_km[i1] > outerStep_km): + Ye[i1]= outerYe + print(tLog.OK+'Y_e Step Function Created.') + return Ye +# Create Ye as a Smooth Function at Inner Area +def smoothInnerYe_for_Athar1995cx(distAll_km, innerStep_km=695.700, outerStep_km=695700\ + , innerYe=0.35, midYe=0.4998, outerYe=0.75): + print(tLog.INFO+'Creating Ye as a Smooth Function at Inner Area.') + Ye= np.zeros(len(distAll_km)) + for i1 in range(len(distAll_km)): + if (distAll_km[i1] <= innerStep_km): + Ye[i1]= ((midYe- innerYe)/ innerStep_km)* distAll_km[i1]+ innerYe + elif (distAll_km[i1] <= outerStep_km and distAll_km[i1] > innerStep_km): + Ye[i1]= midYe + elif (distAll_km[i1] > outerStep_km): + Ye[i1]= outerYe + print(tLog.OK+'Y_e Smooth Function at Inner Area Created.') + return Ye +# ============================ + +# ============================ +# Create Shockwave Added Baryon Density On PreSN Data +def create_shockedNb_g_cm3(dataPreSN_PATH, tPostBounce_s): + print(tLog.INFO+'Creating Shockwave Added Baryon Density On PreSN Data.') + dataSN_t0s= np.loadtxt(dataPreSN_PATH) + r_cm= dataSN_t0s[:, 0] + r_km= r_cm* 1e-5; lenData= len(r_km) + + if tPostBounce_s < 1: + print(tLog.INFO+'Shockwave Added Baryon Density On PreSN Data is not created.') + print(tLog.ERROR+'Consider tPostBounce_s >= 1s.') + exit(tLog.EXIT) + elif tPostBounce_s >= 1: + density_shock_g_cm3 = np.zeros((lenData)) + r_s_km = (-4.6e+3)+ ((11.3e+3)* tPostBounce_s)+ (0.5* (0.2e+3)* tPostBounce_s**2) + + for it1 in range(lenData): + + if r_km[it1] <= r_s_km: + f_x_function = np.exp(0.28- (0.69* np.log(r_s_km))* ((np.arcsin(1- (r_km[it1]/ r_s_km)))**(1.1)) ) + + density_shock_g_cm3[it1]= 10* f_x_function* dataSN_t0s[it1, 2] + elif r_km[it1] > r_s_km: + density_shock_g_cm3[it1]= dataSN_t0s[it1, 2] + print(tLog.OK+'Shockwave Added Baryon Density On PreSN Data Created.') + return density_shock_g_cm3 +# ============================ + +if __name__ == '__main__': + # Current Directory + currDir = osCommand.getcwd() + "/" + COLLECTIVE_NU_OSC_DIR= str(Path(currDir).parent)+ "/" + + # Initialize parser + parser = argparse.ArgumentParser() + + # Adding optional argument + parser.add_argument("-t", "--Time", help= "Time after bounce in seconds", type = float, required=True) + parser.add_argument("-p", "--PathSN", help= "Path to PreSN Data", required=True) + + # Read arguments from command line + args = parser.parse_args() + + # Check for --Time + if args.Time: + tPostBounce_s= float(args.Time) + print(tLog.INFO+'Time after bounce is', tPostBounce_s, 'seconds.') + else: + print(tLog.ERROR+'Time after bounce is not provided.') + exit(tLog.EXIT) + + # Check for --PathSN + if args.PathSN: + dataPreSN_PATH= args.PathSN + print(tLog.INFO+'Path of PreSN Data is', dataPreSN_PATH) + else: + print(tLog.ERROR+'Path of PreSN Data is not provided.') + exit(tLog.EXIT) + + outputPATH= COLLECTIVE_NU_OSC_DIR+"/backgroundProfiles/SN_distCM_Ye_Nb_g_cm3_t"+str(tPostBounce_s)+"s.dat" + + dataSN_t0s= np.loadtxt(dataPreSN_PATH) + r_cm= dataSN_t0s[:, 0] + r_km= dataSN_t0s[:, 0]* 1e-5 + + # Create Shockwave Added Baryon Density On PreSN Data + density_shock_g_cm3= create_shockedNb_g_cm3(dataPreSN_PATH, tPostBounce_s) + + # Create Ye as a Step Function + Ye= stepYe_for_Athar1995cx(r_km) + + # Create Output Data + dataOut= np.zeros([len(r_cm),3]) + for it1 in range(len(r_km)): + dataOut[it1, 0]= r_cm[it1] + dataOut[it1, 1]= Ye[it1] + dataOut[it1, 2]= density_shock_g_cm3[it1] + np.savetxt(outputPATH, dataOut, fmt='%1.6E') + + print(tLog.OK+'Output Data Created at '+outputPATH+'.') + + + diff --git a/scripts/plotGraphsAfterSimulation.py b/scripts/plotGraphsAfterSimulation.py new file mode 100644 index 0000000..8cc0216 --- /dev/null +++ b/scripts/plotGraphsAfterSimulation.py @@ -0,0 +1,79 @@ +""" +-------------------------------------------------------------------------------- +Copyright (c) 2019, Taygun Bulmus +All rights reserved. +See the LICENSE file for license information. +-------------------------------------------------------------------------------- +----- +ID : plotGraphsAfterSimulation.py +Author : Taygun Bulmus +E-Mail : bulmust@gmail.com +----- +""" +# ============================ +import os as osCommand # File operations +# Current Directory +currDir = osCommand.getcwd() + "/" +from pathlib import Path +# collNuPy main directory +COLLECTIVE_NU_OSC_DIR= str(Path(currDir).parent)+ "/" +from sys import path +# Insert Modules Path +path.insert(0, COLLECTIVE_NU_OSC_DIR+'modules') +# Insert plotgraphs.py path +import plotgraphs as plot + +# ============================ +# Colors For Printing +class tLog: + INFO = '\033[94m'+'[INFO] '+'\033[0m' + OK = '\033[92m'+'[OK] '+'\033[0m' + WARNING = '\033[93m'+'[WARNING] '+'\033[0m' + ERROR = '\033[91m'+'[ERROR] '+'\033[0m' + EXIT = '\033[91m' + '[EXIT] plotGraphsAfterSimulation.py' + '\033[0m' +# ============================ + +if __name__ == "__main__": + RESULTS_SIMULATION_DIR= COLLECTIVE_NU_OSC_DIR+ 'results/'+ input('Enter the simulation name: ')+'/' + NPZ_FILE_PATH= RESULTS_SIMULATION_DIR+ 'data/rhoFlavAll.npz' + + # Try to reach the file + print(tLog.INFO+'Try to reach the file:', NPZ_FILE_PATH) + if Path(NPZ_FILE_PATH).is_file(): + print(tLog.OK+'File Found.') + else: + print(tLog.ERROR+'File Not Found.') + exit(tLog.EXIT) + print('') + print(tLog.INFO+'Choose plot(s): Available plots are following') + print(tLog.INFO+'1 = distDiag => Diagonal elements of density matrix to distance.') + print(tLog.INFO+'2 = energyDiag => Diagonal elements of density matrix to energy.') + print(tLog.INFO+'3 = distHamiltDiag => Diagonal elements of Hamiltonian to distance.') + print('') + print(tLog.INFO+'Example: To plot distDiag and distHamiltDiag, enter 13') + plotCheck= input('Enter the plot numbers: ') + print(tLog.INFO+'The plotting is started.') + if '1' in plotCheck: + print(tLog.INFO+'Plotting distDiag.') + plot.plotMain(RESULTS_SIMULATION_DIR).distDiag() + print(tLog.OK+'Plotting distDiag is done.') + if '2' in plotCheck: + print(tLog.INFO+'Plotting energyDiag.') + plot.plotMain(RESULTS_SIMULATION_DIR).energyDiag() + print(tLog.OK+'Plotting energyDiag is done.') + if '3' in plotCheck: + print(tLog.INFO+'Plotting distHamiltDiag.') + plot.plotMain(RESULTS_SIMULATION_DIR).distHamiltDiag() + print(tLog.OK+'Plotting distHamiltDiag is done.') + print(tLog.OK+'The plotting is done.') + + #! ============================ + #! WORKS ONLY FOR UNIX SYSTEMS + # Clear Unused Files and Folders + osCommand.chdir(COLLECTIVE_NU_OSC_DIR) + # Find *.pyc files and remove them + osCommand.system('find . -name "*.pyc" -exec rm -f {} \;') + # Find all __pycache__ folders and remove them + osCommand.system('find . | grep -E "(__pycache__|.pyc|.pyo$)" | xargs rm -rf') + print(tLog.OK+"Unused Files and Folders Removed.") + #! ============================ \ No newline at end of file diff --git a/scripts/random_ri_magPolDecDist/run_random_ri_magPolDecDist.py b/scripts/random_ri_magPolDecDist/run_random_ri_magPolDecDist.py new file mode 100644 index 0000000..50fd84b --- /dev/null +++ b/scripts/random_ri_magPolDecDist/run_random_ri_magPolDecDist.py @@ -0,0 +1,137 @@ +""" +-------------------------------------------------------------------------------- +Copyright (c) 2019, Taygun Bulmus +All rights reserved. +See the LICENSE file for license information. +-------------------------------------------------------------------------------- +----- +ID : plotGraphsAfterSimulation.py +Author : Taygun Bulmus +E-Mail : bulmust@gmail.com +----- +""" +import numpy as np +import random +import os as osCommand # File operations +import argparse + +# Current Directory +currDir = osCommand.getcwd() + "/" +from pathlib import Path +# collNuPy main directory +COLLECTIVE_NU_OSC_DIR= str(Path(str(Path(currDir).parent)+ "/").parent) + '/' + +# ============================ +# Colors For Printing +class tLog: + INFO = '\033[94m'+'[INFO] '+'\033[0m' + OK = '\033[92m'+'[OK] '+'\033[0m' + WARNING = '\033[93m'+'[WARNING] '+'\033[0m' + ERROR = '\033[91m'+'[ERROR] '+'\033[0m' + EXIT = '\033[91m' + '[EXIT] run_random_ri_magPolDecDist.py' + '\033[0m' +# ============================ + +if __name__ =='__main__': + # ============================ + # Initialize parser + parser = argparse.ArgumentParser() + # Adding optional argument + parser.add_argument("-i", "--Ei_MeV", help= "Initial Energy", type= float) + parser.add_argument("-f", "--Ef_MeV", help= "Final Energy", type= float) + parser.add_argument("-m", "--NumberOfMod", help= "Number of Mod", type= int) + # Read arguments from command line + args = parser.parse_args() + # Check for --Ei_MeV + if args.Ei_MeV: + Ei_MeV= float(args.Ei_MeV) + print(tLog.INFO+'Initial Energy is ', Ei_MeV, ' MeV.') + else: + print(tLog.INFO+'Initial Energy is not given. Using default value of 1.0 MeV.') + Ei_MeV= 1.0 + # Check for --Ef_MeV + if args.Ef_MeV: + Ef_MeV= float(args.Ef_MeV) + print(tLog.INFO+'Final Energy is ', Ef_MeV, ' MeV.') + else: + print(tLog.INFO+'Final Energy is not given. Using default value of 50.0 MeV.') + Ef_MeV= 1.0 + # Check for --NumberOfMod + if args.NumberOfMod: + Emod= int(args.NumberOfMod) + print(tLog.INFO+'Number of Energy Mod is ', Emod, '.') + else: + print(tLog.INFO+'Number of Energy Mod is not given. Using default value of 500.') + Emod= 500 + # ============================ + + # Energies + energy= np.zeros(Emod); energy[0]= Ei_MeV + + # Initial Distance (Different For Every Energy) + ri_km= np.zeros(Emod) + + # magneticField_polynomialDecayDistance (Different For Every Energy) + magneticField_polynomialDecayDistance= np.zeros(Emod) + + # Starting Random Points + ri_randStart_km= 49.95 + magneticField_polynomialDecayDistance_randStart= 49.95 + + # Assign First Values + ri_km[0]= ri_randStart_km+ 0.1* random.random() + magneticField_polynomialDecayDistance[0]= magneticField_polynomialDecayDistance_randStart+ 0.1* random.random() + + # Create Random ri_km and magneticField_polynomialDecayDistance + for i1 in range(1, int(Emod)): + # Create Energy + energy[i1]= float(Ei_MeV) + ((float(Ef_MeV)- float(Ei_MeV))/ (Emod- 1))* i1 + + # Assign Random ri_km and magneticField_polynomialDecayDistance + ri_km[i1]= ri_randStart_km+ 0.1* random.random() + magneticField_polynomialDecayDistance[i1]= 49.95+ 0.1* random.random() + + # Start Doing The Simulation + for it_do in range(Emod): + # ============================ + # Modify New ri_km and magneticField_polynomialDecayDistance in physicalParameters.dat + file= open(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters.dat', 'r') + fileWrite= open(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters2.dat', 'w') + for line in file.readlines(): + # Energy + if (line.startswith('numberOf_energyMode')): + line= 'numberOf_energyMode=1 # Valid for only integer\n' + fileWrite.write(line) + elif (line.startswith('energy_initial')): + line= 'energy_initial='+ str(energy[it_do])+ ' # [MeV]\n' + fileWrite.write(line) + elif (line.startswith('energy_final')): + line= 'energy_final='+ str(energy[it_do])+ ' # [MeV]\n' + fileWrite.write(line) + # ri_km + elif (line.startswith('distance_initial')): + line= 'distance_initial='+ str(ri_km[it_do])+ ' # [km]\n' + fileWrite.write(line) + # magneticField_polynomialDecayDistance + elif (line.startswith('magneticField_polynomialDecayDistance')): + line= 'magneticField_polynomialDecayDistance='+ str(magneticField_polynomialDecayDistance[it_do])+ ' # [km] Valid for (use_defaultMagneticProfile=1, magneticField_profile=2)\n' + fileWrite.write(line) + else: + fileWrite.write(line) + file.close() + fileWrite.close() + # Remove and Rename + osCommand.remove(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters.dat') + osCommand.rename(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters2.dat', COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters.dat') + # ============================ + osCommand.chdir(COLLECTIVE_NU_OSC_DIR) + #! ============================ + #! WORKS ONLY FOR UNIX SYSTEMS + osCommand.system('python3 '+COLLECTIVE_NU_OSC_DIR+'collNuPy.py') + #! ============================ + exit() + + + + + + diff --git a/scripts/run_random_ri_magPolDecDist.py b/scripts/run_random_ri_magPolDecDist.py new file mode 100644 index 0000000..3ba24b7 --- /dev/null +++ b/scripts/run_random_ri_magPolDecDist.py @@ -0,0 +1,103 @@ +""" +-------------------------------------------------------------------------------- +Copyright (c) 2019, Taygun Bulmus +All rights reserved. +See the LICENSE file for license information. +-------------------------------------------------------------------------------- +----- +ID : plotGraphsAfterSimulation.py +Author : Taygun Bulmus +E-Mail : bulmust@gmail.com +----- +""" +import numpy as np +import random +import os as osCommand # File operations +# Current Directory +currDir = osCommand.getcwd() + "/" +from pathlib import Path +# collNuPy main directory +COLLECTIVE_NU_OSC_DIR= str(Path(currDir).parent)+ "/" + +# ============================ +# Colors For Printing +class tLog: + INFO = '\033[94m'+'[INFO] '+'\033[0m' + OK = '\033[92m'+'[OK] '+'\033[0m' + WARNING = '\033[93m'+'[WARNING] '+'\033[0m' + ERROR = '\033[91m'+'[ERROR] '+'\033[0m' + EXIT = '\033[91m' + '[EXIT] run_random_ri_magPolDecDist.py' + '\033[0m' +# ============================ + +if __name__ =='__main__': + # Energies + Ei_MeV= 1.0; Ef_MeV= 2.0; Emod=2 + energy= np.zeros(Emod); energy[0]= Ei_MeV + + # Initial Distance (Different For Every Energy) + ri_km= np.zeros(Emod) + + # magneticField_polynomialDecayDistance (Different For Every Energy) + magneticField_polynomialDecayDistance= np.zeros(Emod) + + # Starting Random Points + ri_randStart_km= 49.95 + magneticField_polynomialDecayDistance_randStart= 49.95 + + # Assign First Values + ri_km[0]= ri_randStart_km+ 0.1* random.random() + magneticField_polynomialDecayDistance[0]= 49.95+ 0.1* random.random() + + # Create Random ri_km and magneticField_polynomialDecayDistance + for i1 in range(1, int(Emod)): + # Create Energy + energy[i1]= float(Ei_MeV) + ((float(Ef_MeV)- float(Ei_MeV))/ (Emod- 1))* i1 + + # Assign Random ri_km and magneticField_polynomialDecayDistance + ri_km[i1]= ri_randStart_km+ 0.1* random.random() + magneticField_polynomialDecayDistance[i1]= 49.95+ 0.1* random.random() + + # Start Doing The Simulation + for it_do in range(Emod): + # ============================ + # Modify New ri_km and magneticField_polynomialDecayDistance in physicalParameters.dat + file= open(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters.dat', 'r') + fileWrite= open(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters2.dat', 'w') + for line in file.readlines(): + # Energy + if (line.startswith('numberOf_energyMode')): + line= 'numberOf_energyMode=1 # Valid for only integer\n' + fileWrite.write(line) + elif (line.startswith('energy_initial')): + line= 'energy_initial='+ str(energy[it_do])+ ' # [MeV]\n' + fileWrite.write(line) + elif (line.startswith('energy_final')): + line= 'energy_final='+ str(energy[it_do])+ ' # [MeV]\n' + fileWrite.write(line) + # ri_km + elif (line.startswith('distance_initial')): + line= 'distance_initial='+ str(ri_km[it_do])+ ' # [km]\n' + fileWrite.write(line) + # magneticField_polynomialDecayDistance + elif (line.startswith('magneticField_polynomialDecayDistance')): + line= 'magneticField_polynomialDecayDistance='+ str(magneticField_polynomialDecayDistance[it_do])+ ' # [km] Valid for (use_defaultMagneticProfile=1, magneticField_profile=2)\n' + fileWrite.write(line) + else: + fileWrite.write(line) + file.close() + fileWrite.close() + # Remove and Rename + osCommand.remove(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters.dat') + osCommand.rename(COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters2.dat', COLLECTIVE_NU_OSC_DIR+ 'parameters/physicalParameters.dat') + # ============================ + osCommand.chdir(COLLECTIVE_NU_OSC_DIR) + #! ============================ + #! WORKS ONLY FOR UNIX SYSTEMS + osCommand.system('python3 '+COLLECTIVE_NU_OSC_DIR+'collNuPy.py') + #! ============================ + + + + + +