#NEO-HAZARDS - PYTHON VERSION

import math
import random
import pandas as pd
import time
import os
import matplotlib.pyplot as plt

data = []

#random.random = lambda: 0.5

user_defined_type = input("Enter either asteroid or comet (or press Enter): ")
user_defined_mass = input("Enter the objects's mass(g) manually (or press Enter to calculate it automatically): ")
user_defined_velocity = input("Enter the object's velocity(km/s) manually (or press Enter): ")
user_defined_density = input("Enter the object's density(kg/m^3) manually (or press Enter): ")

#subroutine INITIALIZE
SUMIR = 0
SUME = 0
SUMM = 0
G = 6.673e-8
SUNMASS = 1.992e+33
SUMBFAT = 0
SUMFFAT = 0
SUMTSFAT = 0
SUMGFAT = 0
SUMFATAL = 0
AU = 1.49598e+13 #cm

PLANET = "Earth"
MASS = 5.973e+27 #grams
RADIUS = 6.378e+08 #cm
ZO = 120

RHOZERO = 0.001295
H = float(580000)
SUNDIST = float(1) * AU
PAREA = 4 * math.pi * RADIUS ** 2
EDENSITY = 3 * MASS / (4 * math.pi * RADIUS ** 3)
GRAV = G * MASS / (RADIUS ** 2)
VESCKM = math.sqrt(2 * G * MASS / RADIUS) / float(100000)
VORBKM = math.sqrt(G * SUNMASS / SUNDIST) / float(100000)

LOGTSTEP = 0 #(a time step of 1 year).
TSTEPS = 100 #(default number of time steps is 100).

def mass():
	global M, L
	#subroutine MASS
	L = LOGTSTEP + 0.2171 * math.log(TSTEPS)
	S = -0.6

	if L < 4:
		S = 0.15 * L - 1.2
	elif L < 0:
		S = -1.2
	elif L < -4:
		S = -2.6667 - 0.36667 * L
	elif L < -5:
		S = -0.8333

	LT = LOGTSTEP

	Y = (10 ** (1.6667 * (LT - 2.2)) + 10 ** (0.8333 * (LT - 1))) * random.random() ** (1 / S)
	#M = 9e+12 * Y / (30 + 1.67 * 0.4343 * math.log(Y)) ** 2

	if user_defined_mass.strip():
		try:
			M = float(user_defined_mass)  
		except ValueError:
			print("Invalid mass value")
			exit(1)
	else:
		M = 9e+12 * Y / (30 + 1.67 * 0.4343 * math.log(Y)) ** 2
	return

def line1640():
	global TYPE, IR, FH2O, RHO, CR, TAU, HVAP
	HVAP = 8.08e+10  # 1640
	C = random.random()  # 1640

	if C < 0.082:  # 1640
		if C < 0.031:  # 1710
			if C < 0.021:  # 1730
				IR = 0.0000035  # 1770
				TYPE = "ast/iron"
				FH2O = 0 #water
				RHO = 7.9 #density (g cm-3)
				CR = 3.7e+09 #Strength (dyn cm-2)
				TAU = 0.3 
				C = random.random()  # 1770
				HVAP = 8.26e+10
				if C > 0.995:  # 1880
					IR = 0.000058  # 1880
				elif C > 0.97:
					IR = 1.5e-07
				elif C > 0.935:
					IR = 7.5e-08
				elif C > 0.9:
					IR = 1.5e-08
				elif C > 0.86:
					IR = 0.000035
				elif C > 0.81:
					IR = 3.5e-08
				elif C > 0.74:
					IR = 3.5e-07
				elif C > 0.67:
					IR = 7.49e-07
				elif C > 0.57:
					IR = 0.000015
				elif C > 0.46:
					IR = 0.0000075
				elif C > 0.29: #1780
					IR = 0.0000015 #1780
			else:
				TAU = 0.3  # 1740
				if random.random() < 0.73:  # 1740
					IR = 3e-08  # 1760
					TYPE = "ast/pall"
					FH2O = 0
					RHO = 4.3 + random.random()
					CR = 3e+09  # 1760
				else:
					IR = 0.000002  # 1750
					TYPE = "ast/meso"
					FH2O = 0
					RHO = 5
					CR = 3e+09
		else:
			IR = 1e-08  # 1720
			TYPE = "ast/ach"
			FH2O = 0
			RHO = 3.2
			CR = 1.7e+09
			TAU = 0.02  # 1720
	else:
		if random.random() < 0.77:  # 1650
			if random.random() < 0.72:  # 1680
				IR = 0.0000005  # 1700
				TYPE = "ast/Och "
				FH2O = 0
				RHO = 3.61
				CR = 7e+08
				TAU = 0.2  # 1700
			else:
				IR = 0.0000004  # 1690
				TYPE = "ast/CMch"
				FH2O = 0.1
				RHO = 2.75
				CR = 5e+07
				TAU = 0.05  # 1690
		else:
			IR = 0.0000003  # 1670
			TYPE = "ast/CIch"
			FH2O = 0.2
			RHO = 2.25
			CR = 1e+07
			TAU = 0.05  # 1670
	return

def kind():
	global TYPE, FLAG
	C = random.random()  # 1560
	FLAG = ""  # 1560
	if user_defined_type.strip():
		try:
			if user_defined_type.lower()=="asteroid":
				TYPE = "asteroid"
				return line1640()
			elif user_defined_type.lower()=="comet":
				TYPE = "comet"
				return
			else:
				print("Invalid Type")
				exit(1)
		except ValueError:
			exit(1)
	else:
		if M > 4.5e+16:  # 1560
			RATIO = 10 ** (5.38 - (0.294 * math.log(M) / 2.303))  # 1600
			if M > 10 ** 18:  # 1610
				RATIO = 1  # 1610
			CFRACT = 1 - (RATIO / (1 + RATIO))  # 1610
			if CFRACT > random.random():  # 1620
				TYPE = "comet" #1630
				return # 1630
			else:
				return line1640()
		else:
			if C < 0.75:  # 1570
				TYPE = "asteroid"
				return line1640()
			else:
				TYPE = "comet" # 1630
				return

def energy():
	global TYPE, VIMP, EIMP, MIR, IRID, MH2OIMP, MTON, E, ELEV, R, IR, FH2O, RHO, CR, HVAP,TAU
	#subroutine ENERGY
	if TYPE == "comet": #1920
		HVAP = 3.2e+10  # 1940
		TAU = 0.02
		if random.random() > 0.5:  # 2130
			IR = 0.0000001
			TYPE = "com/long"
			if random.random() > 0.5:  # 2180
				VIMP = (46 - 27 * random.random() ** 2) * 100000.0  # 2200

			else:
				VIMP = (46 + 27 * random.random() ** 2) * 100000.0  # 2190
			FH2O = 0.75
		RHO = 1.1
			CR = 1e+08  # 2210

		else:
			IR = 0.0000002  # 2140
			TYPE = "com/per "
			if random.random() > 0.5:
				VIMP = (24 - 11 * random.random() ** 2) * 100000.0  # 2160
			else:
				VIMP = (24 + 11 * random.random() ** 2) * 100000.0  # 2150
			FH2O = 0.5
			RHO = 1.4
			CR = 1e+08  # 2170

	else:
		VIMP = 100000.0 * (1.01 * VESCKM + VORBKM * (1 - math.cos(random.random()) ** 2))  # 1930

	if user_defined_velocity.strip():
		try:
			VIMP = float(user_defined_velocity)*10**5
		except ValueError:
			print("invalid velocity")
			exit(1)
	else:
		pass
	
	MIR = M * IR  # 2320
	EIMP = M * VIMP ** 2 / 2
	IRID = MIR * 100000.0 / PAREA
	MH2OIMP = FH2O * M

	MTON = EIMP / (4.2e+22)  # 2330
	Q = 2 * random.random() - 1
	ELEV = 28.6479 * (1.570796 - math.atan(Q / math.sqrt(1 - Q * Q)))

	if user_defined_density.strip():
		try:
			RHO = float(user_defined_density)/1000
		except ValueError:
			print("invalid density")
			exit(1)
	else:
		pass
	R = (3 * M / (4 * math.pi * RHO)) ** 0.333  # 2340
	if EIMP > 5e+34:
		E = "e"
	else:
		E = " "

	return

def blowoff():
	global VI, FLAG, MH2OKEPT,  MIR, IRID, VIMP
	#subroutine BLOWOFF
	VI = VIMP / 100000.0 #2380
	if TYPE[:3] == "ast": #2450
		if math.log(M)/2.303 > (20*math.log(VIMP)/2.303-110): #2490
			FB = 0.1 * VI - 1.75
		else:
			FB = -7.63 + 0.276 * math.log(EIMP) / 2.303 #2500
	else:
		if math.log(M)/2.303 > (27.5 * math.log(VIMP)/2.303 - 160): #2460
			FB = 0.13333 * VI - 2.3
		else:
			FB = -5.38 + .2 * math.log(EIMP) / 2.303 #2470

	FB = max(0.0, min(FB, 1.0))

	MIR *= (1 - FB) #2550
	IRID *= (1 - FB) #2560
	MH2OKEPT = MH2OIMP * (1 - FB)

	FLAG = "B" if FB == 1 else "b" if FB > 0 else " "

	return

MM_list = []
V_list = []
Z_list = []

def fragment():
	global F, W, CR, HVAP, Z, RHOATM, TAU, VIMP, ZMIN, V, L, T, MMAX, E, MM ,MO, MM_list, V_list, Z_list

	F = 0
	MO = M
	MM = M #
	L = 0
	T = 0
	FC = 1
	EL = ELEV / 57.2958
	if RHOZERO == 0:
		return

	V = VIMP
	W = 0
	Z = 100000 * ZO
	ZMIN = Z
	VZ = V * math.sin(EL)
	VX = V * math.cos(EL)
	OLDM = 0
	#print(MM)
	MM_list.clear()
	V_list.clear()
	Z_list.clear()
	MM_list.append(MM)
	V_list.append(V)
	Z_list.append(Z)

	for ZSTEP in range(1, 2001):
		R = (3 * MM / (4 * math.pi * RHO)) ** 0.333
		AREA = math.pi * R * R
		CRS = CR * (4.4 / MM) ** 0.0833
		RHOATM = RHOZERO * math.exp(-Z / H)

		DYNP = 0.5 * RHOATM * V * V
		DT = 290000 * FC / V
		if Z < 300000:
			DT = 50000 * FC / V

		DV = -DYNP * AREA * DT / MM
		DVZ = DV * math.sin(EL) * DT / MM + GRAV * DT
		DVX = DV * math.cos(EL)
		LAMBRHO = 1e-08 * (1.37 + 1.8 * math.sqrt(0.00001 * V))
		if V > 2200000:
			LAMBRHO = 1.4e-07
		D = 3 * R * RHOATM / (4 * LAMBRHO)
		GAMMA = 0.6 / math.sqrt(D) + 0.36 / D
		if GAMMA > 1:
			GAMMA = 1
		DM = -0.413 * GAMMA * AREA * RHOATM * DT * V ** 3 / HVAP
		RAT = (DM / MM + DVX / VX + abs(DVZ / VZ))
		if RAT < -0.3:
			FC = -0.1 * FC / RAT
			continue

		DZVEL = -(VZ + 0.5 * DVZ) * DT
		DZGRAV = -GRAV * DT * DT / 2
		DZGEOM = (VX * DT) ** 2 / (RADIUS + Z)
		MM += DM
		DZ = DZVEL + DZGRAV + DZGEOM
		VZ = -DZ / DT
		V = math.sqrt((V + DV) ** 2 - GRAV * DZ)
		MM_list.append(MM)
		V_list.append(V)

		if VZ > V:
			VZ = V
		VX = math.sqrt(V ** 2 - VZ ** 2)
		if Z + DZ < Z:
			ZMIN = Z + DZ
		if ZMIN < 0:
			ZMIN = 0
		if VX == 0:
			VX = 0.000001
		if Z + DZ < 0:
			DT = -DT * Z / DZ #2820
			DZ = -Z
		Z += DZ
		Z_list.append(Z)
		L += V * DT
		T += DT
		EL = math.atan(VZ / VX)
		E = EL * 57.2958
		if Z < 0:
			Z = 0
		if Z != 0:
			
			F = -TAU * DM * V * V / (8 * math.pi * DT * Z * Z) #2850
			MAG = -26.7 - 1.0857 * math.log(F / 1300000.0)

			if MAG < OLDM:
				MMAX = MAG
			OLDM = MAG
			W = 0.5 * MM * V * V

		if MM < 0.0001 * MO: #2880
			MM = 0
			W = 0
			F = "b"
			return

		if V < 300000: #2900
			RRN = random.random()
			if RRN < 0.706:
				F = "od"
			else:
				F = "ld"
			return

		if DYNP > CRS: #2940
			F = "f" #2950
			VDISP = 1.9 * V * math.sqrt(RHOATM / RHO)

			R += VDISP * DT
			AREA = math.pi * R * R
			CRS = 1.2 * CR * (4.4 / MM) ** 0.0833
			if W > 0.5 * MO * VIMP * VIMP:
				continue
			if Z == 0:
				F = "i"
				return
			else:
				return

		if Z < 0:
			F = "i"
			return
		if Z > 100000 * ZO:
			VESC = math.sqrt(2 * G * MASS / (RADIUS + ZO))
			if V < VESC:
				F = "o"
			else:
				F = "e" #should be e
			W = 0
			return
		continue
	return


def crater():
	global RFRAGSTR, RAFF, AAFF, RCRATER, MTON
	RCRATER = 0
	RAFF = 0
	AAFF = 0
	DCRATER = 0
	MEJECT = 0
	MTON = W / 4.2E+22

	RFRAGSTR = H * math.sqrt(RHOZERO / RHO)
	if R > RFRAGSTR or F == "i":
		RCRATER = 4630.0*(1000 * MTON) ** 0.3 #crater radius (diameter)
		RAFF = 6 * RCRATER
		AAFF = math.pi * RAFF * RAFF
	else:
		return

	DCRATER = 2800 * (1000 * MTON) ** 0.3 #crater diameter
	MEJECT = 3.1416 * RCRATER * RCRATER * 2.6 * DCRATER
	MDUST = 0.01 * MEJECT
	return

def hazard():  
    global F,W,Z,R,RFRAGSTR,MTON,RHOATM, FATAL
    #subroutine HAZARD
    FATAL = 0
    LS = -8
    BFATAL = 0
    FFATAL = 0
    TSFATAL = 0
    GFATAL = 0
    LL = random.random()

    LO = -9 + 8 * LL ** 4
    SIGMAOCN = 10 ** LO

    DISTKM = 0
    RUNUPHT = 0
    FLOODPEN = 0
    KK = random.random()
    R4PSIKM = 0
    R1PSIKM = 0

    if F == "e" or F == "b" or F == "o":
    	return
    C = math.log(random.random()) / -0.22
    LS = 10 ** (0.1737 * math.log(C))
    SIGMALND = 0.2 * 10 ** LS  # population density
    YMT = W / 4.185e+22
    MNOX = 1.3e-13 * W
    XNOX = 1000000.0 * MNOX / (RHOZERO * H * PAREA)  # Nox

    ZKM = Z / 100000.0
    YMT3 = YMT ** 0.333
    R4PSIKM = 2.09 * ZKM - 0.45 * ZKM * ZKM / YMT3 + 5.08 * YMT3

    R1PSIKM = 10 * YMT3 * (1.7 * ZKM ** 0.5 - 1.36 * ZKM ** 2 + 1.158)

    if R1PSIKM < 0:
        R1PSIKM = 0
    if R4PSIKM < 0:
        R4PSIKM = 0

    if F == "i" or (F == "f" and R > RFRAGSTR): #goto3330
        RSLANT = 1180000.0*math.sqrt(MTON) #3340
        if KK>0.706: #goto3430	
            F = "li "
            AAFF = math.pi*R4PSIKM**2
            BFATAL = AAFF*SIGMALND
            #goto3580
            GFATAL = 0.3 * SIGMAOCN * R1PSIKM ** 2  # 3580
            if KK > 0.706:
                GFATAL = 0.3 * SIGMALND * R1PSIKM ** 2
            if FFATAL < 0.5:
                FFATAL = 0
            if BFATAL < 0.5:
                BFATAL = 0
            if TSFATAL < 0.5:
                TSFATAL = 0
            if GFATAL < 0.5:
                GFATAL = 0
            FATAL = BFATAL + FFATAL + TSFATAL + GFATAL
            return
        else:
            F = "oi"
            AAFF = math.pi*R4PSIKM**2
            BFATAL = 0.7*AAFF*SIGMAOCN
		
            DISTKM = 3657-5750000.0*math.sqrt(4.04496e-07 - 0.0000004*random.random())
            DISTCM = 100000.0 * DISTKM

            if MTON > 125:
                WVAMPCM = 101500.0*math.sqrt(math.sqrt(MTON))/DISTKM
            else:
                WVAMPCM = 30400.0* math.sqrt (MTON)/DISTKM  

            RUNUPHT = 30*WVAMPCM
            if RUNUPHT>150:
                FLOODPEN=10*(RUNUPHT-150)**(4/3)
            else: 
                FLOODPEN=0
       
            AFLOOD = DISTKM*FLOODPEN/100000.0
            TSFATAL=AFLOOD*SIGMALND

            FATAL=BFATAL+FFATAL+TSFATAL+GFATAL #3640
            return

    else: #goto3440
        if F == "f" and R<RFRAGSTR and Z>0:
            #goto3450
            ZKM = Z / 100000.0  # 3450
            AAFF = math.pi * R4PSIKM ** 2
            if KK > 0.706: #3470
                BFATAL = AAFF * SIGMALND
            if KK >0.706:
                F = "lf"
            if KK < 0.706: #3490
                BFATAL = AAFF * SIGMAOCN
            if KK < 0.706: #3490
                F = "of"
            if KK < 0.706: #3490
                LS = LO

            TFIRE = 0.038 * ((1000 * MTON) ** 0.44)*(RHOATM / RHOZERO) ** 0.36

            RSLANT = 0
            RAFF = 0
            AAFF = 0
            FFATAL = 0

            FLUEXPL = 4.4e+21 * MTON/(math.pi*Z*Z)
            if FLUEXPL > 1e+09:
                RSLANT = ZKM*math.sqrt(FLUEXPL/1e+09)
                RAFF = 0.7*math.sqrt(RSLANT**2-ZKM**2)

                AAFF=math.pi*RAFF**2
                FFATAL=0.5*AAFF*SIGMALND        
                if KK<0.706:
                    FFATAL=0.5*AAFF*SIGMAOCN

            #goto3580
            GFATAL = 0.3 * SIGMAOCN * R1PSIKM ** 2  # 3580
            if KK > 0.706:
                GFATAL = 0.3 * SIGMALND * R1PSIKM ** 2
            if FFATAL < 0.5:
                FFATAL = 0
            if BFATAL < 0.5:
                BFATAL = 0
            if TSFATAL < 0.5:
                TSFATAL = 0
            if GFATAL < 0.5:
                GFATAL = 0
            FATAL = BFATAL + FFATAL + TSFATAL + GFATAL
            return

#main program
def monte_carlo():   #monte_carlo(t_steps)

	simulations = []

	for event in range(1, 745):
		mass() #defines G.V M
		kind() #defines TYPE, FLAG and TYPE, IR, FH2O, RHO, CR, TAU, HVAP
energy()  #energyV defines G.V VIMP, EIMP, MIR, IRID, MH2OIMP, MTON, E, ELEV, R and              redefines TYPE
		blowoff()     #defines G.V VI, FLAG, MH2OKEPT
		fragment()
		crater()
		hazard()

		simulations.append({
			"Event no.": event,
            			"Type": TYPE,  
            			"Mass (kg)": M/1000,
            			"M(H2O)": MH2OIMP/1000,
            	"Velocity (km/s)": VIMP/10**5,
            	"Energy (J)": EIMP/(10*7),
            	"Iridium (ug/cm^2)": IRID,
            	"F": F,
            	"Elev (deg)": ELEV,
            	"Zmin (km)": ZMIN/10**5,
            	"Vfinal (km/s)": V/10**5,
            	"Mmax": MMAX,
			"FATAL": FATAL,
			"Density (kg/m^3)": RHO*1000,
			"HVAP (j/kg)":HVAP/1000,
			"Radius (m)":R/100,
			"Crater (m)":RCRATER/100,
			"RFRAGSTR": RFRAGSTR/100,
			"MASS LEFT (%)": 100*MM/MO
       		 })

		
		#print(V_list)
		perM_lost = [] # percentage of mass lost
		velocity = []
		height = []
		perM_lost.clear()
		velocity.clear()
		height.clear()
		#print(velocity)

		for i in range(len(MM_list)):
			perM_lost.append((MM_list[0]-MM_list[i])*100/MM_list[0])
			velocity.append((V_list[i]/100000))
			height.append(Z_list[i]/100000)

		#print(velocity)
		data.append({
			"event": event,
			"perM_lost": perM_lost,
			"velocity": velocity,
			"height": height
			})

	df1 = pd.DataFrame(simulations)
	script_dir = os.path.dirname(os.path.abspath(__file__))
	timestamp = time.strftime("%Y%m%d_%H%M%S")  # Format: 20241216_142305
	filename1 = f"hazards_{timestamp}.xlsx"
	file_path = os.path.join(script_dir, filename1)
	df1.to_excel(file_path, index=False)

	df2 = pd.DataFrame(data)
	filename2 = f"graphical_data_{timestamp}.xlsx"
	file_path = os.path.join(script_dir, filename2)
	df2.to_excel(file_path, index=False)

	print("Results of Simulation results saved to disk")

monte_carlo()