pip install numpy matplotlib

# calculate_schwarzschild.py
import numpy as np

# Constants
G = 6.674e-11  # Gravitational constant (m^3 kg^-1 s^-2)
c = 2.998e8    # Speed of light (m/s)

def calculate_schwarzschild_radius(mass_kg):
    """
    Calculates the Schwarzschild radius for a given mass.

    Args:
        mass_kg (float): Mass of the object in kilograms.

    Returns:
        float: Schwarzschild radius in meters.
    """
    if mass_kg <= 0:
        raise ValueError("Mass must be positive.")
    
    Rs = (2 * G * mass_kg) / (c**2)
    return Rs

if __name__ == "__main__":
    # Example 1: Earth's mass
    M_earth = 5.972e24  # Mass of Earth in kg
    Rs_earth = calculate_schwarzschild_radius(M_earth)
    print(f"Schwarzschild Radius for Earth ({M_earth:.2e} kg): {Rs_earth:.4e} meters")
    print(f"This is approximately {Rs_earth * 1000:.2f} millimeters.")

    print("-" * 30)

    # Example 2: Sun's mass
    M_sun = 1.989e30  # Mass of Sun in kg
    Rs_sun = calculate_schwarzschild_radius(M_sun)
    print(f"Schwarzschild Radius for Sun ({M_sun:.2e} kg): {Rs_sun:.4e} meters")
    print(f"This is approximately {Rs_sun / 1000:.2f} kilometers.")

    print("-" * 30)

    # Example 3: A supermassive black hole (e.g., Sagittarius A*)
    # Sagittarius A* (center of Milky Way) is about 4.3 million solar masses
    M_sag_A_star = 4.3e6 * M_sun 
    Rs_sag_A_star = calculate_schwarzschild_radius(M_sag_A_star)
    print(f"Schwarzschild Radius for Sagittarius A* ({M_sag_A_star:.2e} kg): {Rs_sag_A_star:.2e} meters")
    print(f"This is approximately {Rs_sag_A_star / (1.5e11):.2f} AU (Astronomical Units).")
    print(f"Note: 1 AU is approx. 150 million km, Earth-Sun distance.")

python calculate_schwarzschild.py

# black_hole_visualization.py
import numpy as np
import matplotlib.pyplot as plt

# --- Constants (from previous script) ---
G = 6.674e-11  # Gravitational constant (m^3 kg^-1 s^-2)
c = 2.998e8    # Speed of light (m/s)

def calculate_schwarzschild_radius(mass_kg):
    return (2 * G * mass_kg) / (c**2)

# --- Simulation Parameters ---
IMAGE_SIZE = 800  # Resolution of the output image (pixels x pixels)
MASS_OF_BH_SOLAR = 10  # Mass of our simulated black hole in solar masses
M_sun = 1.989e30       # Mass of Sun in kg

# Calculate Rs for our simulated black hole
mass_bh_kg = MASS_OF_BH_SOLAR * M_sun
Rs_actual = calculate_schwarzschild_radius(mass_bh_kg) # In meters

# Scale factor for our visualization:
# How many pixels does 1 meter represent? Let Rs_actual take up a certain fraction of the image.
# We want Rs_actual to be visible, so let's normalize our coordinates.
# Let the width of our simulation space be, say, 10 * Rs_actual.
sim_width = 10 * Rs_actual # meters
pixel_scale = IMAGE_SIZE / sim_width # pixels per meter

# Normalized Rs in simulation units (pixels for drawing)
# Rs in terms of image pixels.
# The apparent horizon (photon sphere) is typically around 1.5 * Rs.
# Everything inside this radius on screen will be black.
# We map the center of the image to (0,0) of the black hole.
normalized_Rs = Rs_actual * pixel_scale
apparent_horizon_radius_pixels = 2.5 * normalized_Rs # A slightly larger visual "shadow"

print(f"Simulating a {MASS_OF_BH_SOLAR} solar mass black hole.")
print(f"Actual Schwarzschild Radius: {Rs_actual:.2e} meters")
print(f"Normalized Schwarzschild Radius (pixels): {normalized_Rs:.2f} pixels")
print(f"Apparent Horizon Radius (pixels): {apparent_horizon_radius_pixels:.2f} pixels")
print("-" * 40)

# --- Create Background Image (a simple grid) ---
def create_grid_background(size, num_lines=20):
    background = np.zeros((size, size, 3), dtype=np.uint8) # RGB image
    line_spacing = size // num_lines
    
    # Draw horizontal lines
    for i in range(num_lines):
        y = i * line_spacing
        background[y:y+2, :, :] = 200 # Light gray lines
    
    # Draw vertical lines
    for i in range(num_lines):
        x = i * line_spacing
        background[:, x:x+2, :] = 200 # Light gray lines
    
    # Make background darker
    background[background == 0] = 50 
    
    return background

background_image = create_grid_background(IMAGE_SIZE)

# --- Initialize Output Image ---
output_image = np.copy(background_image) # Start with background

# --- Gravitational Lensing Calculation ---
# Create coordinate grids
# X and Y go from -sim_width/2 to sim_width/2 in meters
x_coords_meters = np.linspace(-sim_width / 2, sim_width / 2, IMAGE_SIZE)
y_coords_meters = np.linspace(-sim_width / 2, sim_width / 2, IMAGE_SIZE)
X_meters, Y_meters = np.meshgrid(x_coords_meters, y_coords_meters)

# Calculate radial distance in meters
R_meters = np.sqrt(X_meters**2 + Y_meters**2)

# Calculate angles (for polar coordinates)
Theta = np.arctan2(Y_meters, X_meters)

# Define the "pull" factor for lensing. This is a visual approximation.
# Objects closer to the black hole are displaced more.
# K is a constant to tune the strength of the lensing effect.
# A value around 0.5 to 2.0 often works well for visual effect.
K_lensing = 1.0 

# Apparent (visual) radius of the black hole's shadow. 
# Inside this, we simply render black.
# This is typically larger than Rs due to strong lensing effects.
# The photon sphere is at 1.5 * Rs, but the "shadow" is larger.
# Let's use the 'apparent_horizon_radius_pixels' already calculated for our screen.
# Convert back to meters for comparison with R_meters
apparent_horizon_radius_meters = apparent_horizon_radius_pixels / pixel_scale

# Iterate over each pixel in the image
for py in range(IMAGE_SIZE):
    for px in range(IMAGE_SIZE):
        # Current pixel's coordinates in meters relative to black hole center
        current_r_meters = R_meters[py, px]
        current_theta = Theta[py, px]

        # If inside the apparent horizon, it's black
        if current_r_meters < apparent_horizon_radius_meters:
            output_image[py, px] = [0, 0, 0] # Black
        else:
            # Calculate source coordinates for lensing effect
            # We want to find (r_source, theta_source) that maps to (current_r_meters, current_theta)
            # The 'pull' is proportional to K * Rs_actual / r_meters
            
            # Simple radial displacement model:
            # Shift the source outwards radially. This causes background objects to spread out near BH.
            # This is a common way to visualize lensing, effectively "stretching" the background.
            # The larger current_r_meters is, the smaller the pull.
            # Avoid division by zero close to BH (already handled by apparent_horizon check)
            
            # This formula ensures that as r_meters gets closer to Rs_actual, the pull becomes stronger.
            # r_prime is the radial distance in the background image
            
            # Note: The visual effect often comes from displacing the *source* radially outwards.
            # The 2.0 multiplier here is a visual fudge factor to make the effect more pronounced.
            # The 'deflection' should cause the light to appear from further away radially.
            deflection = K_lensing * (Rs_actual**2) / current_r_meters # A simple inverse-r deflection magnitude
            
            r_source_meters = current_r_meters + deflection # Source is further out radially
            
            # Clamp r_source_meters to reasonable bounds to prevent out-of-bounds access
            # This is particularly important for areas very close to the black hole where deflection can be extreme
            max_r_source = sim_width * 0.7 # Limit source to 70% of sim_width to avoid extreme edge cases
            if r_source_meters > max_r_source:
                r_source_meters = max_r_source
            
            # Convert source polar coordinates back to Cartesian (in meters)
            sx_meters = r_source_meters * np.cos(current_theta)
            sy_meters = r_source_meters * np.sin(current_theta)
            
            # Map source meters coordinates back to background image pixel coordinates
            # Center of image (px=IMAGE_SIZE/2, py=IMAGE_SIZE/2) maps to (0,0) meters
            source_px = int(sx_meters * pixel_scale + IMAGE_SIZE / 2)
            source_py = int(sy_meters * pixel_scale + IMAGE_SIZE / 2)
            
            # Ensure source pixel coordinates are within bounds of the background image
            source_px = np.clip(source_px, 0, IMAGE_SIZE - 1)
            source_py = np.clip(source_py, 0, IMAGE_SIZE - 1)
            
            output_image[py, px] = background_image[source_py, source_px]


# --- Draw the actual Schwarzschild Radius (Event Horizon) ---
# It's a perfect circle at normalized_Rs
center_x, center_y = IMAGE_SIZE // 2, IMAGE_SIZE // 2
radius_pixels = int(normalized_Rs)

# Create a circle for the event horizon
theta_circle = np.linspace(0, 2*np.pi, 300)
x_circle = center_x + radius_pixels * np.cos(theta_circle)
y_circle = center_y + radius_pixels * np.sin(theta_circle)

# Draw the circle onto the output image
# Use a bright color (e.g., red or green) for visibility
for i in range(len(x_circle)):
    cx, cy = int(x_circle[i]), int(y_circle[i])
    if 0 <= cx < IMAGE_SIZE and 0 <= cy < IMAGE_SIZE:
        output_image[cy, cx] = [255, 0, 0] # Red line for event horizon

# A slightly thicker line
for i in range(len(x_circle)):
    cx, cy = int(x_circle[i]), int(y_circle[i])
    for dx in [-1, 0, 1]:
        for dy in [-1, 0, 1]:
            nx, ny = cx + dx, cy + dy
            if 0 <= nx < IMAGE_SIZE and 0 <= ny < IMAGE_SIZE:
                output_image[ny, nx] = [255, 0, 0] # Red line for event horizon

# --- Display the image ---
plt.figure(figsize=(8, 8))
plt.imshow(output_image)
plt.title(f"Black Hole Simulation ({MASS_OF_BH_SOLAR} Solar Masses)")
plt.axis('off') # Hide axes ticks
plt.show()

python black_hole_visualization.py