ASS.py Simplified Cosmos Explorer 0.2A

import numpy as np
import matplotlib.pyplot as plt
import pubchempy as pcp
import requests
import json
import time
import re
import warnings
import socket
import struct
import threading
import traceback
import logging
import argparse

# Suppress Astroquery's InsecureRequestWarning when connecting to some services
from requests.packages.urllib3.exceptions import InsecureRequestWarning
warnings.simplefilter('ignore', InsecureRequestWarning)

# --- Library Imports for Astrophysics and Data ---
# This block handles imports for the client-side astrophysics application.
from pystac_client import Client
from datetime import datetime
from astroquery.gaia import Gaia
from astroquery.cadc import Cadc
from astroquery.utils.tap import TapPlus
from astroquery.ned import Ned
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table
import nmrglue as ng # Used by the server for NMR data handling

# --- Core Scientific & Mathematical Constants ---
# This section consolidates constants from all programs.
# Physics Constants (from astro.py)
CURIE_TO_BQ = 3.7e10
LN2 = np.log(2)
STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8
PLANCK_CONSTANT = 6.62607015e-34
BOLTZMANN_CONSTANT = 1.380649e-23
C = 299792458
HBAR = 1.054571817e-34
LAMBDA_C = HBAR / (np.pi * C)

# Networking & Operation Codes (from n-dim.py and NMRBBS.py)
OPERATION_INTERPOLATE = 0
OPERATION_DIFFERENTIATE = 1
OPERATION_CALCULATE_GRADIENT_1D = 2
OPERATION_HYPERBOLIC_INTERCEPT_HANDLER = 3
OPERATION_INTEGRATE = 4
OPERATION_INTEGRATE_ND = 5
OPERATION_WORKFLOW = 6
SERVER_HOST = '127.0.0.1'
SERVER_PORT = 12345

# Server Message Codes (from NMRBBS.py)
MSG_TYPE_INTERPOLATION_RESULT = 5
MSG_TYPE_DIFFERENTIATION_RESULT = 12
MSG_TYPE_CLIENT_DISCONNECT = 254
MSG_TYPE_ID_ASSIGNMENT = 255

# --- Shared Networking & Data Handling Suite ---
# These helper functions are crucial for communication between the client and server.
# They are used by both the `Server` class and the client's `execute_workflow` function.

def _sendall_data(sock, data_array):
    """Helper to send a numpy array's bytes."""
    data_bytes = data_array.astype(np.float32).tobytes()
    sock.sendall(struct.pack('!I', len(data_bytes)))
    sock.sendall(data_bytes)

def _recvall(sock, n):
    """Helper function to reliably receive exactly 'n' bytes from a socket."""
    data = b''
    while len(data) < n:
        packet = sock.recv(n - len(data))
        if not packet:
            return None
        data += packet
    return data

def _recvall_data(sock):
    """Helper to receive a numpy array's bytes or an error message."""
    data_len_bytes = _recvall(sock, 4)
    if data_len_bytes is None:
        raise ConnectionResetError("Incomplete data (length) received from server.")
    data_len = struct.unpack('!I', data_len_bytes)[0]
    data_bytes = _recvall(sock, data_len)
    if data_bytes is None:
        raise ConnectionResetError("Incomplete data (content) received from server.")
    try:
        error_message = data_bytes.decode('utf-8')
        if "Error" in error_message or "Server internal error" in error_message:
            raise RuntimeError(f"Server returned an error: {error_message}")
    except UnicodeDecodeError:
        pass
    return np.array(struct.unpack(f'!{data_len // 4}f', data_bytes), dtype=np.float32)

# --- Server-Side Mathematical & Data Services Suite ---
# This class contains the logic that runs when the program is started in server mode.
# It handles requests for complex mathematical operations and data broadcasting.

logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s')

class Server:
    def __init__(self, host: str, port: int):
        self.host = host
        self.port = port
        self.clients: list[dict] = []
        self.clients_lock = threading.Lock()
        self.client_data: dict[socket.socket, dict] = {}
        self.client_data_lock = threading.Lock()
        self.next_client_id = 0

    def pseudo_interpolate_arcsecant_1d_triple(self, fx: np.ndarray, fy: np.ndarray, fz: np.ndarray, x_interp: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        """
        Performs a custom, non-linear interpolation on 1D data with y and z datasets.
        This function uses a simple sine-based easing function to approximate the non-linear
        interpolation intended in the original code, avoiding common mathematical errors.
        """
        if not (len(fx) == len(fy) == len(fz)) or len(fx) < 2:
            raise ValueError("X, Y, and Z data must have equal length and at least two points.")
        
        interp_y = np.zeros_like(x_interp, dtype=float)
        interp_z = np.zeros_like(x_interp, dtype=float)

        for i, x in enumerate(x_interp):
            if x <= fx[0]:
                idx1, idx2 = 0, 1
            elif x >= fx[-1]:
                idx1, idx2 = len(fx) - 2, len(fx) - 1
            else:
                idx1 = np.searchsorted(fx, x, side='right') - 1
                idx2 = idx1 + 1

            x1, x2 = fx[idx1], fx[idx2]
            y1, y2 = fy[idx1], fy[idx2]
            z1, z2 = fz[idx1], fz[idx2]

            if (x2 - x1) == 0:
                interp_y[i] = y1
                interp_z[i] = z1
                continue

            t_linear = (x - x1) / (x2 - x1)
            scaled_t = 0.5 * (1 + np.sin(np.pi * (t_linear - 0.5)))
            scaled_t = np.clip(scaled_t, 0, 1)

            interp_y[i] = y1 + (y2 - y1) * scaled_t
            interp_z[i] = z1 + (z2 - z1) * scaled_t

        return interp_y, interp_z

    def pseudo_interpolate_arcsecant_nd_triple(self, all_fy_data: list[np.ndarray], all_fz_data: list[np.ndarray], all_fx_data: list[np.ndarray], x_interp: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        """
        Performs the non-linear interpolation for n-dimensional data by calling
        the 1D interpolation function for each dimension.
        """
        all_interp_y = []
        all_interp_z = []
        num_dimensions = len(all_fy_data)

        if not (len(all_fx_data) == num_dimensions == len(all_fz_data)):
            raise ValueError("The number of x, y, and z data arrays must match across dimensions.")

        for i, (fx, fy, fz) in enumerate(zip(all_fx_data, all_fy_data, all_fz_data)):
            interp_y, interp_z = self.pseudo_interpolate_arcsecant_1d_triple(fx, fy, fz, x_interp)
            all_interp_y.extend(interp_y)
            all_interp_z.extend(interp_z)

        return np.array(all_interp_y), np.array(all_interp_z)

    def numerical_derivative_1d(self, y_values: np.ndarray, x_values: np.ndarray) -> np.ndarray:
        """
        Calculates the numerical derivative of a 1D array using a robust gradient method.
        """
        if len(y_values) != len(x_values) or len(y_values) < 2:
            raise ValueError("Y and X data must have equal length and at least two points for derivative calculation.")
        return np.gradient(y_values, x_values)

    def differentiate_arcsecant_nd_triple(self, all_fy_data: list[np.ndarray], all_fz_data: list[np.ndarray], all_fx_data: list[np.ndarray], x_eval: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
        """
        Calculates the numerical derivative for n-dimensional data by differentiating
        the interpolated data for each dimension.
        """
        num_dimensions = len(all_fy_data)
        if not (len(all_fx_data) == num_dimensions == len(all_fz_data)):
            raise ValueError("The number of x, y, and z data arrays must match.")
        if len(x_eval) < 2:
            raise ValueError("At least two evaluation points are needed for differentiation.")

        all_derivatives_y = []
        all_derivatives_z = []

        for i, (fx, fy, fz) in enumerate(zip(all_fx_data, all_fy_data, all_fz_data)):
            interpolated_y, interpolated_z = self.pseudo_interpolate_arcsecant_1d_triple(fx, fy, fz, x_eval)
            derivatives_y = self.numerical_derivative_1d(interpolated_y, x_eval)
            derivatives_z = self.numerical_derivative_1d(interpolated_z, x_eval)
            all_derivatives_y.extend(derivatives_y)
            all_derivatives_z.extend(derivatives_z)

        return np.array(all_derivatives_y), np.array(all_derivatives_z)

    def handle_client(self, client_socket: socket.socket, client_address: tuple):
        """Main loop for handling client connections."""
        logging.info(f"Connection from {client_address}")
        try:
            # We only expect a single workflow request in this simplified example
            op_code_byte = _recvall(client_socket, 1)
            if not op_code_byte: return
            op_code = struct.unpack('!B', op_code_byte)[0]
            
            if op_code == OPERATION_WORKFLOW:
                self.handle_workflow_request(client_socket)
        except Exception as e:
            logging.error(f"Error handling client {client_address}: {e}", exc_info=True)
        finally:
            client_socket.close()

    def handle_workflow_request(self, client_socket):
        """Processes a workflow request from a client."""
        logging.info("Processing workflow request...")
        data_len_bytes = _recvall(client_socket, 4)
        if not data_len_bytes: return
        data_len = struct.unpack('!I', data_len_bytes)[0]
        data_bytes = _recvall(client_socket, data_len)
        if not data_bytes: return
        
        workflow_steps = json.loads(data_bytes.decode('utf-8'))
        
        # This is where the magic happens: the server processes the client's request.
        # It's a simplified version for demonstration.
        try:
            # For this integrated version, we assume the workflow only contains the
            # interpolation and differentiation steps, as defined in the client code.
            interpolated_data = None
            for step in workflow_steps:
                if step["operation_type"] == "INTERPOLATE":
                    input_data = step["input_data"]
                    fx = np.array(input_data["fx_data"][0])
                    fy = np.array(input_data["fy_data"][0])
                    fz = np.array(input_data["fz_data"][0])
                    x_interp = np.array(step["parameters"]["x_interp_points"])
                    interp_y, interp_z = self.pseudo_interpolate_arcsecant_nd_triple([fy], [fz], [fx], x_interp)
                    interpolated_data = np.concatenate([interp_y, interp_z])
                elif step["operation_type"] == "DIFFERENTIATE":
                    # The client's workflow references the previous step's output
                    input_data = step["input_data"]
                    if input_data.get("source_id") == "interpolated_data" and interpolated_data is not None:
                        # Split the interpolated data back into y and z for differentiation
                        mid_point = len(interpolated_data) // 2
                        interp_y = interpolated_data[:mid_point]
                        interp_z = interpolated_data[mid_point:]
                        x_eval = np.array(step["parameters"]["x_eval_points"])
                        
                        # Note: The server's differentiate function handles the rest of the logic
                        # by re-interpolating and differentiating.
                        deriv_y, deriv_z = self.differentiate_arcsecant_nd_triple([interp_y], [interp_z], [x_eval, x_eval], x_eval)
                        
                        # Send back the result
                        result_array = np.concatenate([deriv_y, deriv_z])
                        _sendall_data(client_socket, result_array)
                        logging.info("Successfully processed and sent back differentiation result.")
                        return

        except Exception as e:
            error_msg = f"Server internal error: {e}"
            logging.error(error_msg, exc_info=True)
            try:
                client_socket.sendall(struct.pack('!I', len(error_msg.encode('utf-8'))))
                client_socket.sendall(error_msg.encode('utf-8'))
            except: pass


    def run(self):
        """Starts the server and listens for connections."""
        server_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
        server_socket.bind((self.host, self.port))
        server_socket.listen(5)
        logging.info(f"Server listening on {self.host}:{self.port}")
        while True:
            try:
                client_socket, client_address = server_socket.accept()
                client_thread = threading.Thread(target=self.handle_client, args=(client_socket, client_address))
                client_thread.daemon = True
                client_thread.start()
            except KeyboardInterrupt:
                logging.info("Server shutting down.")
                break
            except Exception as e:
                logging.error(f"Server accept loop error: {e}")
        server_socket.close()

# --- Client-Side Application & Physics Exploration Suite ---
# This section contains the user-facing application logic. When run in client mode,
# it explores various physics concepts and makes requests to the server for calculations.

def execute_workflow(workflow_steps):
    """
    Executes a sequence of operations on the server as a single workflow.
    This is the bridge between the client's physics exploration and the server's math services.
    """
    client_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
    try:
        client_socket.connect((SERVER_HOST, SERVER_PORT))
        print(f"Connected to server at {SERVER_HOST}:{SERVER_PORT}")
        client_socket.sendall(struct.pack('!B', OPERATION_WORKFLOW))
        workflow_json = json.dumps(workflow_steps)
        workflow_bytes = workflow_json.encode('utf-8')
        client_socket.sendall(struct.pack('!I', len(workflow_bytes)))
        client_socket.sendall(workflow_bytes)
        print("Workflow sent to server. Waiting for result...")
        result = _recvall_data(client_socket)
        return result
    except Exception as e:
        print(f"Client error during workflow execution: {e}")
        return None
    finally:
        client_socket.close()
        print("Connection closed.")

def _encode_output_with_math(x_data, y_data):
    """
    This function sends data to the server for processing.
    It takes our data points, smooths them out with a line, and then asks the server
    to calculate the slope at different points, which helps us understand how the data is changing.
    """
    print("\n--- Applying advanced math to analyze our data via the server ---")
    
    x_eval_points = np.linspace(min(x_data), max(x_data), 10, dtype=np.float32)

    workflow = [
        {
            "operation_type": "INTERPOLATE",
            "input_data": {
                "type": "direct",
                "fx_data": [x_data.tolist()],
                "fy_data": [y_data.tolist()],
                "fz_data": [[]]
            },
            "parameters": {
                "x_interp_points": x_eval_points.tolist()
            },
            "output_id": "interpolated_data"
        },
        {
            "operation_type": "DIFFERENTIATE",
            "input_data": {
                "type": "reference",
                "source_id": "interpolated_data"
            },
            "parameters": {
                "x_eval_points": x_eval_points.tolist()
            }
        }
    ]

    result = execute_workflow(workflow)
    if result is not None:
        print("Encoded Output (The 'rate of change' of our data from the server):")
        print(result)
        return result
    return None

def analyze_fundamental_physics():
    """Explores fundamental physics concepts and uses the math server for analysis."""
    print("\n--- Exploring Foundational Physics Concepts ---")
    
    mass_kg = 1.0
    mass_energy = mass_kg * C**2
    print(f"A single kilogram of matter contains an immense amount of energy: {mass_energy:.4e} Joules.")

    T_surface = 5778
    solar_radiation = STEFAN_BOLTZMANN_CONSTANT * T_surface**4
    print(f"The Sun's surface at {T_surface} Kelvin emits energy at a rate of: {solar_radiation:.4e} Watts per square meter.")
    
    x_data = np.array([1, 2, 3, 4, 5], dtype=np.float32)
    y_data = np.array([10, 12, 15, 19, 25], dtype=np.float32)
    
    _encode_output_with_math(x_data, y_data)


def explore_cosmic_stability():
    """Analyzes the stability of cosmic objects."""
    print("\n--- Analyzing the Stability of Cosmic Objects ---")

    def _check_stability(object_type, stability_data):
        """Checks if a cosmic object is stable based on its data."""
        print(f"Checking the stability of a '{object_type}'...")
        for point, value in stability_data.items():
            if value < 0:
                print(f"  Alert: Instability detected! A negative value was found at {point}.")
            else:
                print(f"  Status: Stable. Everything looks good at {point}.")

    neutron_star_data = {'pressure_point_A': 10, 'pressure_point_B': -5}
    _check_stability("Neutron Star", neutron_star_data)

    white_dwarf_data = {'pressure_point_A': 5, 'pressure_point_B': 8}
    _check_stability("White Dwarf", white_dwarf_data)


def find_celestial_objects():
    """Queries various astronomical databases."""
    print("\n--- Searching International Astronomical Databases ---")
    
    ra, dec = 180.0, -60.0
    target_name = "M104"
    start_date = "2023-01-01"
    end_date = "2023-12-31"
    
    print(f"Searching for information around coordinates RA:{ra}, DEC:{dec}...")
    
    try:
        coord = SkyCoord(ra=ra, dec=dec, unit=(u.degree, u.degree))
        gaia_data = Gaia.query_object_async(coord, radius=u.Quantity(1.0, u.deg))
        if gaia_data and isinstance(gaia_data, Table):
            print(f"  ✅ Found {len(gaia_data)} stars and objects in the Gaia catalog.")
        else:
            print("  ❌ No Gaia data found or query failed.")
    except Exception as e:
        print(f"  ❌ Gaia search failed: {e}")

    try:
        cadc = Cadc()
        cadc_result = cadc.query_object(target_name, collection='CFHT')
        if cadc_result and len(cadc_result) > 0:
            print(f"  ✅ Found {len(cadc_result)} observations of '{target_name}' in the CADC database.")
        else:
            print(f"  ❌ No CADC data found for '{target_name}'.")
    except Exception as e:
        print(f"  ❌ CADC search failed: {e}")

    try:
        ned_result = Ned.query_object(target_name)
        if ned_result and len(ned_result) > 0:
            print(f"  ✅ Found {len(ned_result)} entries for '{target_name}' in the NED database.")
        else:
            print(f"  ❌ No NED data found for '{target_name}'.")
    except Exception as e:
        print(f"  ❌ NED search failed: {e}")

    try:
        catalog = Client.open("https://earth-search.aws.element84.com/v1")
        search = catalog.search(
            intersects={"type": "Point", "coordinates": [ra, dec]},
            datetime=f"{start_date}T00:00:00Z/{end_date}T23:59:59Z",
        )
        item_collection = search.item_collection()
        print(f"  ✅ Found {len(item_collection)} catalog items in the STAC search.")
    except Exception as e:
        print(f"  ❌ STAC search failed: {e}")
        
    try:
        isro_data_url = "https://isro.gov.in/launcher"
        requests.get(isro_data_url)
        print("  ✅ Successfully connected to the ISRO data endpoint.")
    except requests.exceptions.RequestException as e:
        print(f"  ❌ ISRO data query failed: {e}")

# --- Main Program Entry Point ---
if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Integrated Client/Server for Astrophysics and Data Analysis.")
    parser.add_argument('--mode', type=str, choices=['server', 'client'], required=True,
                        help='Specify the program mode: "server" to run the mathematical service, or "client" to run the astrophysics explorer.')
    args = parser.parse_args()

    if args.mode == 'server':
        print("--- Starting Mathematical Service Server ---")
        server = Server(SERVER_HOST, SERVER_PORT)
        server.run()
    elif args.mode == 'client':
        print("--- Starting Simplified Astrophysics Explorer (Client) ---")
        print("This program runs a full suite of astronomical analyses and database queries for you.")
        print("It will attempt to connect to the math service server if one is running.")
        
        analyze_fundamental_physics()
        explore_cosmic_stability()
        find_celestial_objects()
        
        print("\n--- All client explorations and searches are complete. ---")

astinc.py Theoretical Astronomical Explorer and Cosmological Incongruency Elucidation Engine 1.1

Please enjoy this theoretical astronomical matter exploration tool, some findings may be hypothetically reasonable.

import numpy as np
import matplotlib.pyplot as plt
import pubchempy as pcp
import requests
import json
import time
import re
import warnings

# Suppress Astroquery's InsecureRequestWarning when connecting to some services
from requests.packages.urllib3.exceptions import InsecureRequestWarning
warnings.simplefilter('ignore', InsecureRequestWarning)

# --- NEW: Import pystac-client for STAC database query ---
from pystac_client import Client
from datetime import datetime

# --- Import Astroquery modules for international databases ---
from astroquery.gaia import Gaia
from astroquery.cadc import Cadc
from astroquery.utils.tap import TapPlus
from astroquery.ned import Ned
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table

# --- Scientific Constants ---
CURIE_TO_BQ = 3.7e10  # Conversion factor: 1 Curie = 3.7 x 10^10 Becquerels
LN2 = np.log(2)  # Natural logarithm of 2, used in half-life calculations

STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8  # Stefan-Boltzmann constant in W/(m^2·K^4)
PLANCK_CONSTANT = 6.62607015e-34  # Planck's constant in J·s
BOLTZMANN_CONSTANT = 1.380649e-23  # Boltzmann constant in J/K
C = 299792458  # Speed of light in m/s
HBAR = 1.054571817e-34  # Reduced Planck constant in J·s
LAMBDA_C = HBAR / (np.pi * C)  # Compton wavelength factor

# --- Relativistic Effects & Quantum-Gravity Bridge ---

def calculate_relativistic_comparison():
    """Calculates and prints relativistic and classical energy comparisons."""
    print("\n--- Relativistic Effects on Matter & Radiation ---")
    print("This function compares different energy models and their effects.")

    mass_kg = 1.0  # kg
    c_speed = 299792458  # m/s

    # Relativistic energy
    def _relativistic_energy(m, c):
        return m * c**2

    # Mass-energy equivalence
    mass_energy = _relativistic_energy(mass_kg, c_speed)
    print(f"Mass-energy equivalence for a {mass_kg} kg object: {mass_energy:.4e} Joules")

    # Blackbody radiation
    def _blackbody_radiation(T):
        return STEFAN_BOLTZMANN_CONSTANT * T**4

    T_surface = 5778  # Temperature of the Sun's surface in Kelvin
    solar_radiation = _blackbody_radiation(T_surface)
    print(f"Blackbody radiation from a surface at {T_surface} K: {solar_radiation:.4e} W/m^2")

def calculate_quantum_gravity_bridge():
    """Performs a calculation for the hypothetical quantum-gravity bridge."""
    print("\n--- Hypothetical Quantum-Gravity Bridge ---")
    print("This section explores a conceptual link between quantum and gravitational phenomena.")

    def _corroborate_gamma_data(gamma_type, data):
        """Mock function to simulate data corroboration."""
        print(f"Corroborating data for {gamma_type}...")
        for source, value in data.items():
            if abs(value - 1.0) > 0.1:
                print(f"  Warning: Data from {source} deviates significantly.")
            else:
                print(f"  Success: Data from {source} is within expected range.")

    test_data = {
        'source_A': 1.01,
        'source_B': 0.98,
        'source_C': 1.15
    }
    _corroborate_gamma_data("quantum_gravity_gamma", test_data)

# --- Equation of State (EOS) & Stability ---

def analyze_eos_stability():
    """Analyzes and prints information on Equation of State (EOS) and stability."""
    print("\n--- Equation of State (EOS) & Stability ---")

    def _analyze_eos_stability_data(eos_type, stability_data):
        """Analyzes a single EOS data set."""
        print(f"Analyzing {eos_type} EOS for stability...")
        for point, value in stability_data.items():
            if value < 0:
                print(f"  Warning: Instability detected at {point} with value {value:.2f}.")
            else:
                print(f"  Stability confirmed at {point} with value {value:.2f}.")

    neutron_star_data = {'pressure_point_A': 10, 'pressure_point_B': -5}
    _analyze_eos_stability_data("neutron_star", neutron_star_data)

    white_dwarf_data = {'pressure_point_A': 5, 'pressure_point_B': 8}
    _analyze_eos_stability_data("white_dwarf", white_dwarf_data)


# --- Consolidated Database Queries ---

def fetch_gaia_data(ra, dec, radius):
    """Fetches data from the Gaia database."""
    print(f"\n--- Querying Gaia for RA:{ra}, DEC:{dec} ---")
    try:
        coord = SkyCoord(ra=ra, dec=dec, unit=(u.degree, u.degree))
        job = Gaia.query_object_async(coord, radius=u.Quantity(radius, u.deg))
        if job and isinstance(job, Table):
            print(f"Found {len(job)} Gaia sources.")
        else:
            print("No Gaia data found or query failed.")
    except Exception as e:
        print(f"Gaia query failed: {e}")

def fetch_cadc_data(target_name):
    """Fetches data from the CADC database."""
    print(f"\n--- Querying CADC for target '{target_name}' ---")
    try:
        cadc = Cadc()
        cadc.query_object(target_name)
        result = cadc.query_object(target_name, collection='CFHT')
        if result and len(result) > 0:
            print(f"Found {len(result)} CADC observations for {target_name}.")
        else:
            print("No CADC data found or query failed.")
    except Exception as e:
        print(f"CADC query failed: {e}")

def fetch_ned_data(target_name):
    """Fetches data from the NED database."""
    print(f"\n--- Querying NED for target '{target_name}' ---")
    try:
        result_table = Ned.query_object(target_name)
        if result_table and len(result_table) > 0:
            print(f"Found {len(result_table)} NED entries for {target_name}.")
        else:
            print("No NED data found or query failed.")
    except Exception as e:
        print(f"NED query failed: {e}")

def fetch_stac_data(ra, dec, start_date, end_date):
    """Fetches STAC catalog data."""
    print("\n--- Querying STAC Catalog ---")
    try:
        catalog = Client.open("https://earth-search.aws.element84.com/v1")
        search = catalog.search(
            intersects={"type": "Point", "coordinates": [ra, dec]},
            datetime=f"{start_date}T00:00:00Z/{end_date}T23:59:59Z",
        )
        item_collection = search.item_collection()
        print(f"Found {len(item_collection)} STAC items.")
    except Exception as e:
        print(f"STAC query failed: {e}")

def fetch_isro_data(target_name):
    """Simulates fetching data from a hypothetical ISRO database."""
    print(f"\n--- Querying Hypothetical ISRO Database for '{target_name}' ---")
    isro_data_url = "https://isro.gov.in/launcher"
    try:
        response = requests.get(isro_data_url)
        response.raise_for_status()
        print(f"Successfully connected to ISRO data endpoint.")
    except requests.exceptions.RequestException as e:
        print(f"ISRO data query failed: {e}")

def query_all_astronomical_databases():
    """Runs a consolidated query against all astronomical databases."""
    print("\n--- Committing Queries Against All Databases ---")
    
    # Example parameters
    ra = 180.0
    dec = -60.0
    target_name = "M104"
    start_date = "2023-01-01"
    end_date = "2023-12-31"
    
    fetch_gaia_data(ra, dec, 1.0)
    fetch_cadc_data(target_name)
    fetch_ned_data(target_name)
    fetch_stac_data(ra, dec, start_date, end_date)
    fetch_isro_data(target_name)

# --- Main Entry Point ---

if __name__ == "__main__":
    print("--- Astrophysics & EOS Stability Explorer (Simplified) ---")
    print("This program runs a full suite of analyses and database queries.")
    
    # Run the consolidated relativistic and quantum physics analysis
    calculate_relativistic_comparison()
    calculate_quantum_gravity_bridge()
    
    # Run the EOS analysis
    analyze_eos_stability()
    
    # Run the consolidated database queries
    query_all_astronomical_databases()

    print("\n--- All analyses and queries complete. ---")

                    

Astronomy Physics Explorer 2.0

Had to roll it back a bit, working under inverse leverages.

import numpy as np
import matplotlib.pyplot as plt
import pubchempy as pcp
import requests
import json
import time
import re
import warnings

# Suppress Astroquery's InsecureRequestWarning when connecting to some services
from requests.packages.urllib3.exceptions import InsecureRequestWarning
warnings.simplefilter('ignore', InsecureRequestWarning)

# --- NEW: Import pystac-client for STAC database query ---
from pystac_client import Client
from datetime import datetime

# --- Import Astroquery modules for international databases ---
from astroquery.gaia import Gaia
from astroquery.cadc import Cadc
from astroquery.utils.tap import TapPlus
from astropy import units as u
from astropy.coordinates import SkyCoord
from astropy.table import Table

# --- Scientific Constants ---
# Nuclear Physics Constants
CURIE_TO_BQ = 3.7e10  # Conversion factor: 1 Curie = 3.7 x 10^10 Becquerels (disintegrations per second)
LN2 = np.log(2)  # Natural logarithm of 2, used in half-life calculations

# Blackbody Radiation Constants
STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8  # Stefan-Boltzmann constant in W/(m^2·K^4)
PLANCK_CONSTANT = 6.62607015e-34  # Planck's constant in J·s
BOLTZMANN_CONSTANT = 1.380649e-23  # Boltzmann constant in J/K
SPEED_OF_LIGHT = 2.99792458e8  # Speed of light in a vacuum in m/s
WIEN_DISPLACEMENT_CONSTANT = 2.897771955e-3  # Wien's displacement law constant in m·K

# General Relativity Constants
GRAVITATIONAL_CONSTANT = 6.67430e-11  # G in m³·kg⁻¹·s⁻²
SPEED_OF_LIGHT_SQ = SPEED_OF_LIGHT**2
SOLAR_MASS_KG = 1.989e30  # Mass of the Sun in kg
PARSEC_M = 3.086e16  # Parsec in meters

# --- Manual Nuclear Data for Gamma Emission ---
MANUAL_GAMMA_DATA = {
    "Cesium-137": {
        "half_life_days": 30.08 * 365.25,
        "gamma_energies_MeV": [0.662],
        "gamma_yields_per_decay": [0.851],
    },
    "Cobalt-60": {
        "half_life_days": 1925.21,
        "gamma_energies_MeV": [1.173, 1.332],
        "gamma_yields_per_decay": [0.999, 0.999],
    },
    "Iodine-131": {
        "half_life_days": 8.02,
        "gamma_energies_MeV": [0.364],
        "gamma_yields_per_decay": [0.817],
    },
    "Potassium-40": {
        "half_life_days": 1.251e9 * 365.25,
        "gamma_energies_MeV": [1.4608],
        "gamma_yields_per_decay": [0.1067],
    },
    "Aluminum-26": {
        "half_life_days": 7.17e5 * 365.25,
        "gamma_energies_MeV": [1.809],
        "gamma_yields_per_decay": [1.0],
    },
    "Iron-60": {
        "half_life_days": 2.6e6 * 365.25,
        "gamma_energies_MeV": [1.173, 1.332],
        "gamma_yields_per_decay": [1.0, 1.0],
    },
}

# --- Astronomical Object Presets ---
# Storing characteristic mass and a typical observation distance
ASTRO_OBJECT_PRESETS = {
    "Neutron Star": {
        "mass_kg": 1.4 * SOLAR_MASS_KG,
        "distance_m": 1.0e5,  # 100 km from center
        "description": "A dense remnant of a massive star.",
    },
    "Stellar Black Hole": {
        "mass_kg": 10 * SOLAR_MASS_KG,
        "distance_m": 1.0e6,  # 1000 km from center
        "description": "A black hole formed from a collapsing star.",
    },
    "Galactic Center (Sgr A*)": {
        "mass_kg": 4.3e6 * SOLAR_MASS_KG,
        "distance_m": 100 * PARSEC_M,  # 100 parsecs away
        "description": "The supermassive black hole at the center of the Milky Way.",
    },
    "Galaxy Cluster": {
        "mass_kg": 1e15 * SOLAR_MASS_KG,
        "distance_m": 1e6 * PARSEC_M,  # 1 Megaparsec into the cluster
        "description": "A massive gravitationally bound collection of galaxies.",
    },
    "Intra-galactic Void": {
        "mass_kg": 1e12 * SOLAR_MASS_KG,  # Mass of a single galaxy as a reference point
        "distance_m": 10e6
        * PARSEC_M,  # 10 Megaparsecs away from the nearest large mass
        "description": "The vast, relatively empty space between galaxies.",
    },
}

# ---------------------------------------------------------------------------
# PubChem Helper Functions
def _get_compound_cid_pubchem(compound_name):
    try:
        compounds = pcp.get_compounds(compound_name, "name")
        return compounds[0].cid if compounds else None
    except Exception:
        return None

def _fetch_half_life_from_pubchem(cid):
    base_url = "https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/"
    url = f"{base_url}cid/{cid}/JSON/"
    try:
        response = requests.get(url, timeout=5)
        response.raise_for_status()
        data = response.json()
        search_sections = [
            "Atomic Mass, Half Life, and Decay",
            "Chemical and Physical Properties",
            "Depositor-Supplied",
        ]
        for record in data.get("Record", {}).get("Section", []):
            if record.get("TOCHeading") in search_sections:
                for prop_section in record.get("Section", []):
                    if "half-life" in prop_section.get("TOCHeading", "").lower():
                        for item in prop_section.get("Information", []):
                            for value_dict in item.get("Value", {}).get(
                                "StringWithMarkup", []
                            ):
                                text = value_dict.get("String")
                                if text:
                                    match = re.match(
                                        r"(\d+\.?\d*(?:e[+-]?\d+)?)\s*(seconds|minutes|hours|days|years|s|min|hr|d|y)",
                                        text.lower(),
                                    )
                                    if match:
                                        return float(match.group(1)), match.group(2)
        return None, None
    except Exception:
        return None, None

def _convert_half_life_to_days(value, unit):
    if unit in ["s", "seconds"]:
        return value / (24 * 3600)
    if unit in ["min", "minutes"]:
        return value / (24 * 60)
    if unit in ["hr", "hours"]:
        return value / 24
    if unit in ["d", "days"]:
        return value
    if unit in ["y", "years"]:
        return value * 365.25
    return None

# ---------------------------------------------------------------------------
# Core Radiation Calculation Functions
def calculate_decay_constant(half_life_days):
    if half_life_days is None or half_life_days <= 0:
        return 0
    half_life_seconds = half_life_days * 24 * 3600
    return LN2 / half_life_seconds

def calculate_activity_at_time(
    initial_activity_Bq, decay_constant, time_elapsed_seconds
):
    if decay_constant == 0:
        return initial_activity_Bq
    return initial_activity_Bq * np.exp(-decay_constant * time_elapsed_seconds)

def calculate_blackbody_spectral_radiance(wavelength_m, temperature_K):
    if temperature_K <= 0 or wavelength_m <= 0:
        return 0.0
    term1 = (2 * PLANCK_CONSTANT * SPEED_OF_LIGHT_SQ) / (wavelength_m**5)
    term2_exp_arg = (PLANCK_CONSTANT * SPEED_OF_LIGHT) / (
        wavelength_m * BOLTZMANN_CONSTANT * temperature_K
    )
    if term2_exp_arg > 700:
        return 0.0
    spectral_radiance = term1 / (np.exp(term2_exp_arg) - 1)
    return spectral_radiance

# ---------------------------------------------------------------------------
# General Relativity Functions
def calculate_schwarzschild_radius(mass_kg):
    """Calculates the Schwarzschild radius for a given mass."""
    if mass_kg <= 0:
        return 0.0
    return (2 * GRAVITATIONAL_CONSTANT * mass_kg) / SPEED_OF_LIGHT_SQ

def calculate_dilation_factor(mass_kg, distance_m):
    """
    Calculates the time dilation factor for an observer at infinity relative to an object in a gravitational field.
    A factor of 1 means no dilation. Values < 1 mean time is dilated (slowed down).
    """
    if mass_kg <= 0 or distance_m <= 0:
        return 1.0  # No dilation for zero mass or distance
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    if distance_m <= schwarzschild_radius:
        return 0.0  # At or inside event horizon, time stops for an outside observer
    return np.sqrt(1 - (schwarzschild_radius / distance_m))

def calculate_relativistic_activity(
    isotope_name, initial_activity_Bq, time_elapsed_days, mass_kg, distance_m
):
    """
    Calculates activity as observed by a distant observer, accounting for gravitational time dilation
    and a hypothetical "time crystal" effect.
    """
    half_life_days = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("half_life_days")
    if half_life_days is None:
        cid = _get_compound_cid_pubchem(isotope_name)
        if cid:
            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)
            if hl_val:
                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)
    if half_life_days is None:
        return None

    time_elapsed_seconds = time_elapsed_days * 24 * 3600

    # Apply gravitational time dilation
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

    # An observer far away sees time running slower for the object, so the decay process appears slowed.
    effective_time_local_seconds = time_elapsed_seconds * dilation_factor

    # --- Hypothetical: curvature-based modifier ---
    curvature_parameter = calculate_gravitational_curvature_parameter(
        mass_kg, distance_m
    )
    time_crystal_modifier = _calculate_time_crystal_half_life_modifier(
        curvature_parameter
    )

    effective_half_life_days = half_life_days * time_crystal_modifier
    decay_constant = calculate_decay_constant(effective_half_life_days)
    current_activity_Bq = calculate_activity_at_time(
        initial_activity_Bq, decay_constant, effective_time_local_seconds
    )

    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get(
        "gamma_yields_per_decay", [0]
    )
    total_gamma_yield_per_decay = sum(gamma_yields)
    total_gamma_disintegrations_per_second = (
        current_activity_Bq * total_gamma_yield_per_decay
    )
    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ
    return gamma_curies

def calculate_relativistic_peak_wavelength(local_temperature_K, mass_kg, distance_m):
    """Calculates the observed peak wavelength of a blackbody, accounting for gravitational redshift."""
    if local_temperature_K <= 0:
        return None
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    observed_temperature_K = local_temperature_K * dilation_factor
    if observed_temperature_K <= 0:
        return np.inf
    peak_wavelength_m = WIEN_DISPLACEMENT_CONSTANT / observed_temperature_K
    return peak_wavelength_m * 1e9  # Convert to nm

# ---------------------------------------------------------------------------
# Quantum and Gravitational Bridge Functions (hypothetical)
def calculate_gravitational_curvature_parameter(mass_kg, distance_m):
    """
    A simplified, unitless parameter representing local gravitational curvature.
    Higher values imply greater curvature.
    """
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    if distance_m <= 0 or distance_m <= schwarzschild_radius:
        return np.inf
    return schwarzschild_radius / distance_m

def quantum_differentiation_hypothetical(curvature_parameter, base_quantum_value):
    """
    Hypothetical function to illustrate 'quantum mechanical differentiation'.
    NOT based on established physics.
    """
    if curvature_parameter >= 1:
        return 0
    differentiated_value = base_quantum_value * (1 - curvature_parameter)
    return differentiated_value

def _calculate_time_crystal_half_life_modifier(curvature_parameter):
    """
    Hypothetical function to model a 'time crystal' effect driven by curvature.
    Returns a multiplicative modifier on the half-life.
    """
    if curvature_parameter >= 1:
        return 0.0
    modifier = np.tanh(1 - curvature_parameter)
    return modifier

# ---------------------------------------------------------------------------
# EOS & Stability (Nebulas + Cosmic Strings)
PROTON_MASS = 1.67262192369e-27
RADIATION_CONSTANT = 4 * STEFAN_BOLTZMANN_CONSTANT / SPEED_OF_LIGHT  # a_rad
K_B = BOLTZMANN_CONSTANT

# ---------- Gas / Nebula EOS ----------
def sound_speed_ideal_gas(T_K, mu=2.33, gamma=5 / 3):
    """
    Isentropic sound speed for an ideal gas:
      c_s = sqrt(gamma * k_B * T / (mu * m_p))
    mu ~ 2.33 for molecular clouds (H2 + He).
    """
    return np.sqrt(gamma * K_B * T_K / (mu * PROTON_MASS))

def pressure_ideal_gas(rho, T_K, mu=2.33):
    return rho * K_B * T_K / (mu * PROTON_MASS)

def pressure_polytrope(rho, K, gamma):
    return K * rho**gamma

def sound_speed_polytrope(rho, K, gamma):
    # c_s^2 = dP/dρ = gamma * K * ρ^(gamma-1)
    return np.sqrt(max(0.0, gamma * K * rho ** (gamma - 1)))

def pressure_radiation(T_K):
    # P_rad = (1/3) a T^4
    return (1.0 / 3.0) * RADIATION_CONSTANT * T_K**4

def mixed_barotropic_pressure(rho, T_K, mu=2.33, gamma=5 / 3):
    # simple additive gas + radiation pressure
    return pressure_ideal_gas(rho, T_K, mu) + pressure_radiation(T_K)

# ---------- Nebula stability ----------
def jeans_length(rho, T_K, mu=2.33, gamma=5 / 3):
    """
    λ_J = c_s * sqrt(pi / (G * ρ))
    """
    cs = sound_speed_ideal_gas(T_K, mu, gamma)
    return cs * np.sqrt(np.pi / (GRAVITATIONAL_CONSTANT * rho))

def jeans_mass(rho, T_K, mu=2.33, gamma=5 / 3):
    """
    M_J ~ (4π/3) * (λ_J/2)^3 * ρ
    """
    lam = jeans_length(rho, T_K, mu, gamma)
    return (4 * np.pi / 3.0) * (lam / 2.0) ** 3 * rho

def bonnor_ebert_mass(T_K, P_ext, mu=2.33):
    """
    Critical Bonnor–Ebert mass for isothermal sphere:
      M_BE ≈ 1.18 * (c_s^4) / (G^(3/2) * P_ext^(1/2))
    """
    cs = sound_speed_ideal_gas(T_K, mu, gamma=1.0)  # isothermal -> gamma=1
    return 1.18 * cs**4 / (GRAVITATIONAL_CONSTANT**1.5 * np.sqrt(P_ext))

# ---------- Cosmic strings (coarse-grained) ----------
def string_wiggliness_params(alpha=0.0, mu_base=1.0):
    """
    Phenomenology:
      U = μ (1+α),      T = μ /(1+α)
    α >= 0 encodes small-scale structure ("wiggles").
    """
    U = mu_base * (1.0 + alpha)
    T = mu_base / (1.0 + alpha)
    return U, T

def string_wave_speeds(alpha=0.0):
    """
    For the model above:
      c_T^2 = T/U = 1/(1+α)^2
    Use c_L ~ 1 (Nambu–Goto) as a simple default.
    """
    c_T = 1.0 / (1.0 + alpha)
    c_L = 1.0  # simple Nambu–Goto proxy
    return c_T, c_L

def string_effective_w(v_rms=0.5):
    """
    Effective EOS parameter for a network (toy model):
      w ≈ -1/3 * (1 - v_rms^2)
    """
    v2 = np.clip(v_rms**2, 0.0, 1.0)
    return -1.0 / 3.0 * (1.0 - v2)

def string_loop_oscillation_time(loop_length_m, alpha=0.0):
    """
    Fundamental oscillation time τ ~ L / (c_T * c).
    """
    c_T, _ = string_wave_speeds(alpha)
    return loop_length_m / (c_T * SPEED_OF_LIGHT)

# ---------- Visualizers ----------
def plot_nebula_stability_grid(
    rho_vals, T_vals, mu=2.33, gamma=5 / 3, mass_scale_Msun=1.0
):
    """
    Heatmap: M_J / (mass_scale) to see where collapse is favored (M_J smaller than cloud mass).
    """
    from matplotlib.colors import LogNorm

    R, T = np.meshgrid(rho_vals, T_vals)
    MJ = np.vectorize(jeans_mass)(R, T, mu, gamma)
    plt.figure(figsize=(9, 6))
    im = plt.pcolormesh(R, T, MJ / (mass_scale_Msun * 1.98847e30), norm=LogNorm())
    plt.colorbar(im, label=r"$M_J/M_\odot$")
    plt.xscale("log")
    plt.yscale("log")
    plt.xlabel(r"Density $\rho$ (kg m$^{-3}$)")
    plt.ylabel("Temperature (K)")
    plt.title("Nebula Stability Map (Jeans Mass in Solar Masses)")
    plt.tight_layout()
    plt.show()

def plot_string_speeds_vs_wiggliness(alpha_vals):
    cT = [string_wave_speeds(a)[0] for a in alpha_vals]
    plt.figure(figsize=(8, 5))
    plt.plot(alpha_vals, cT)
    plt.xlabel("Wiggliness α")
    plt.ylabel("Transverse wave speed c_T (units of c)")
    plt.title("Cosmic String Wave Speed vs Wiggliness")
    plt.grid(True)
    plt.tight_layout()
    plt.show()

# ---------------------------------------------------------------------------
# Analysis Functions
def generate_normal_distribution_series(mean, std_dev, num_points=200, spread_factor=3):
    start = mean - spread_factor * std_dev
    end = mean + spread_factor * std_dev
    if start < 0:
        start = 0
    series_values = np.linspace(start, end, num_points)
    pdf_values = (1 / (std_dev * np.sqrt(2 * np.pi))) * np.exp(
        -0.5 * ((series_values - mean) / std_dev) ** 2
    )
    return series_values, pdf_values

def get_gamma_curies_single_point(isotope_name, initial_activity_Bq, time_elapsed_days):
    half_life_days = None
    if isotope_name in MANUAL_GAMMA_DATA:
        half_life_days = MANUAL_GAMMA_DATA[isotope_name]["half_life_days"]
    else:
        cid = _get_compound_cid_pubchem(isotope_name)
        if cid:
            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)
            if hl_val:
                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)
    if half_life_days is None:
        return None
    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get(
        "gamma_yields_per_decay", [0]
    )
    total_gamma_yield_per_decay = sum(gamma_yields)
    decay_constant = calculate_decay_constant(half_life_days)
    time_elapsed_seconds = time_elapsed_days * 24 * 3600
    current_activity_Bq = calculate_activity_at_time(
        initial_activity_Bq, decay_constant, time_elapsed_seconds
    )
    total_gamma_disintegrations_per_second = (
        current_activity_Bq * total_gamma_yield_per_decay
    )
    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ
    return gamma_curies

def analyze_gamma_curies_comparison(
    isotope_name,
    initial_activity_Bq,
    mean_time_days,
    std_dev_time_days,
    mass_kg,
    distance_m,
    num_points=200,
    object_name="Massive Object",
):
    """
    Compares gamma curies decay in flat spacetime vs. near a massive object.
    """
    print(
        f"\n--- Comparing Gamma Curies: Flat Spacetime vs. Relativistic Environment ---"
    )
    print(f"  Isotope: {isotope_name}, Initial Activity: {initial_activity_Bq:.2e} Bq")
    print(
        f"  Relativistic Environment: {object_name} (Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m)"
    )
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")
    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")
    if dilation_factor < 0.99:
        print(
            f"  (A dilation factor less than 1 means time runs slower in the gravitational field.)"
        )
    else:
        print(f"  (Dilation factor is near 1, indicating a weak gravitational effect.)")

    time_series_days, _ = generate_normal_distribution_series(
        mean_time_days, std_dev_time_days, num_points=num_points
    )

    # Flat Spacetime Calculation
    flat_curies_series = [
        get_gamma_curies_single_point(isotope_name, initial_activity_Bq, t_days)
        for t_days in time_series_days
    ]

    # Relativistic Calculation
    relativistic_curies_series = [
        calculate_relativistic_activity(
            isotope_name, initial_activity_Bq, t_days, mass_kg, distance_m
        )
        for t_days in time_series_days
    ]

    if all(c is None for c in flat_curies_series):
        print(
            f"  Calculation failed. Please ensure '{isotope_name}' has half-life data."
        )
        return

    # Plotting
    plt.figure(figsize=(12, 7))
    plt.plot(
        time_series_days, flat_curies_series, label="Flat Spacetime (Local Observer)"
    )
    plt.plot(
        time_series_days,
        relativistic_curies_series,
        label=f"Near {object_name} (Distant Observer)",
        linestyle="--",
    )
    plt.title(f"Gamma Curies Decay: {isotope_name} - Flat vs. {object_name}")
    plt.xlabel("Time (days)")
    plt.ylabel("Gamma Curies (Ci)")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

def simulate_blackbody_wien_comparison(
    temperatures_K, mass_kg, distance_m, object_name="Massive Object"
):
    """
    Compares blackbody peak wavelength in flat spacetime vs. near a massive object.
    """
    print(f"\n--- Comparing Blackbody Peak Wavelength: Flat vs. Relativistic ---")
    print(f"  Analyzing temperatures: {temperatures_K} K")
    print(
        f"  Relativistic Environment: {object_name} (Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m)"
    )
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")
    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")

    temp_list_sorted = sorted(temperatures_K)
    flat_peak_wl_list_nm = [
        (WIEN_DISPLACEMENT_CONSTANT / T) * 1e9 for T in temp_list_sorted if T > 0
    ]
    relativistic_peak_wl_list_nm = [
        calculate_relativistic_peak_wavelength(T, mass_kg, distance_m)
        for T in temp_list_sorted
        if T > 0
    ]

    # Plotting
    plt.figure(figsize=(8, 6))
    plt.plot(
        temp_list_sorted,
        flat_peak_wl_list_nm,
        "o-",
        label="Flat Spacetime (Local Observer)",
    )
    plt.plot(
        temp_list_sorted,
        relativistic_peak_wl_list_nm,
        "r--",
        label=f"Near {object_name} (Distant Observer)",
    )
    plt.title(f"Wien's Law: Peak Wavelength vs. Temperature (Flat vs. {object_name})")
    plt.xlabel("Temperature (K)")
    plt.ylabel("Peak Wavelength (nm)")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

def explore_quantum_gravitational_differentiation(
    mass_kg, distance_m_series, object_name="Massive Object"
):
    """
    New function to explore the hypothetical relationship between gravity and a quantum state.
    """
    print("\n--- Exploring a Hypothetical Gravitational-Quantum Bridge ---")
    print(f"  Analyzing Quantum State near a {object_name} with Mass={mass_kg:.2e} kg")
    curvature_series = [
        calculate_gravitational_curvature_parameter(mass_kg, d)
        for d in distance_m_series
    ]
    base_quantum_value = 1.0
    differentiated_value_series = [
        quantum_differentiation_hypothetical(c, base_quantum_value)
        for c in curvature_series
    ]

    plt.figure(figsize=(10, 6))
    plt.plot(distance_m_series, differentiated_value_series, "b-")
    plt.axvline(
        x=calculate_schwarzschild_radius(mass_kg),
        color="r",
        linestyle="--",
        label="Schwarzschild Radius",
    )
    plt.title(f'Hypothetical "Quantum Differentiation" vs. Distance from {object_name}')
    plt.xlabel(f"Distance from {object_name} (m)")
    plt.ylabel("Hypothetical Differentiated Quantum Value (Unitless)")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

# --- Database Fetching Functions ---
DARTS_TAP_URL = 'https://darts.isas.jaxa.jp/tap/'
darts_tap = TapPlus(url=DARTS_TAP_URL)
NADC_API_BASE_URL = 'http://nadc.china-vo.org/api/query'
STAC_API_URL = 'https://stac.astrogeology.usgs.gov/stac' # Example public STAC endpoint

def fetch_gaia_data(source_id: int) -> Table:
    print(f"🔎 Querying ESA Gaia database for Source ID: {source_id}")
    try:
        query = f"SELECT TOP 1 * FROM gaiadr3.gaia_source WHERE source_id={source_id}"
        job = Gaia.launch_job(query)
        results = job.get_results()
        print("✅ Gaia Data Found:"); print(results)
        return results
    except Exception as e:
        print(f"❌ Error fetching Gaia data: {e}")
        return None

def fetch_neossat_data(coordinates: SkyCoord, radius: u.Quantity = 1 * u.arcsec) -> Table:
    print(f"🔎 Querying CADC for NEOSSAT data around {coordinates.to_string('hmsdms')}")
    try:
        result_table = Cadc.query_region(coordinates, radius=radius, collection='neossat')
        print("✅ NEOSSAT Data Found:"); print(result_table)
        return result_table
    except Exception as e:
        print(f"❌ Error querying CADC for NEOSSAT data: {e}")
        return None

def fetch_darts_data(target_name: str) -> Table:
    print(f"🔎 Querying JAXA DARTS for: {target_name}")
    try:
        adql_query = f"SELECT TOP 10 * FROM TAP_SCHEMA.tables WHERE table_name LIKE '%{target_name}%'"
        job = darts_tap.launch_job_async(adql_query)
        results = job.get_results()
        print("✅ DARTS Data Found:"); print(results)
        return results
    except Exception as e:
        print(f"❌ Error querying DARTS via ADQL: {e}")
        return None

def fetch_nadc_data(search_query: str) -> dict:
    print(f"🔎 Querying NADC for: {search_query}")
    params = {'q': search_query, 'format': 'json'}
    try:
        response = requests.get(NADC_API_BASE_URL, params=params, timeout=10)
        response.raise_for_status()
        data = response.json()
        print("✅ NADC Data Found:"); print(json.dumps(data, indent=2))
        return data
    except requests.exceptions.RequestException as e:
        print(f"❌ An error occurred while fetching data from NADC: {e}")
    except json.JSONDecodeError as e:
        print(f"❌ Failed to decode JSON response from NADC: {e}")
    return None

def fetch_stac_data(ra: float, dec: float, start_date: str, end_date: str, collections: list = None):
    print(f"🔎 Querying STAC catalog at {STAC_API_URL} for coordinates ({ra}, {dec}) between {start_date} and {end_date}")
    try:
        catalog = Client.open(STAC_API_URL)
        search = catalog.search(
            datetime=f"{start_date}T00:00:00Z/{end_date}T23:59:59Z",
            bbox=[ra - 0.05, dec - 0.05, ra + 0.05, dec + 0.05],
            collections=collections
        )
        items = list(search.get_all_items())
        if items:
            print(f"✅ Found {len(items)} STAC Items.")
            for item in items:
                print(f"   - Item ID: {item.id}, Collection: {item.collection_id}, Date: {item.datetime.strftime('%Y-%m-%d')}")
            return items
        else:
            print("❗ No STAC items found for this query.")
            return None
    except Exception as e:
        print(f"❌ An error occurred while fetching STAC data: {e}")
        return None

def fetch_isro_data(target_name: str):
    """
    A placeholder function to demonstrate how to integrate a new database.
    This would be customized based on a public ISRO API, such as for AstroSat data.
    """
    print(f"🔎 Attempting to query ISRO database for: {target_name}")
    try:
        # This is a conceptual API endpoint. A real one would be provided by ISRO's ISSDC.
        # For example, an API for AstroSat data might be structured differently.
        api_url = f"https://api.isro.gov.in/astrosat/data?q={target_name}"
        response = requests.get(api_url, timeout=10)
        response.raise_for_status()
        data = response.json()
        
        # Check for a successful data payload.
        if data and data.get("results"):
            print("✅ ISRO data found:")
            print(json.dumps(data, indent=2))
        else:
            print("❗ No results found from ISRO database.")

    except requests.exceptions.RequestException as e:
        print(f"❌ An error occurred while fetching data from ISRO API: {e}")
        print("  Note: This is a placeholder for a public API that may not be available or may require specific query parameters.")

# --- UI & Main Program Execution ---
def run_relativistic_comparison_menu():
    print("\n--- Relativistic Phenomena Simulation ---")
    print("Select an astronomical environment:")
    for i, (name, data) in enumerate(ASTRO_OBJECT_PRESETS.items()):
        print(f"  {i+1}) {name}: {data['description']}")
    print(f"  {len(ASTRO_OBJECT_PRESETS)+1}) Custom Object")
    try:
        env_choice = int(input("Enter your choice: "))
        if 1 <= env_choice <= len(ASTRO_OBJECT_PRESETS):
            preset_name = list(ASTRO_OBJECT_PRESETS.keys())[env_choice - 1]
            astro_object = ASTRO_OBJECT_PRESETS[preset_name]
            object_mass, object_distance, object_name = astro_object["mass_kg"], astro_object["distance_m"], preset_name
        elif env_choice == len(ASTRO_OBJECT_PRESETS) + 1:
            object_mass = float(input("\nEnter mass of custom object in kg: "))
            object_distance = float(input("Enter distance from the object in meters: "))
            object_name = "Custom Object"
        else:
            print("Invalid choice."); return
        sub_choice = input("Select: (N)uclear Decay or (B)lackbody Radiation? ").strip().lower()
        if sub_choice == 'n':
            isotope_name = input("Enter isotope (e.g., Cesium-137): ").strip()
            initial_activity = float(input("Enter initial activity in Bq: "))
            mean_time = float(input("Enter mean measurement time (days): "))
            std_dev_time = float(input("Enter standard deviation for time (days): "))
            analyze_gamma_curies_comparison(isotope_name, initial_activity, mean_time, std_dev_time, object_mass, object_distance, object_name=object_name)
        elif sub_choice == 'b':
            temp_input_str = input("Enter temperatures (K) separated by commas (e.g., 300, 1000, 5000): ")
            temperatures_list = [float(t.strip()) for t in temp_input_str.split(',')]
            simulate_blackbody_wien_comparison(temperatures_list, object_mass, object_distance, object_name=object_name)
        else:
            print("Invalid choice."); return
    except ValueError:
        print("Invalid input.")

def run_quantum_gravity_menu():
    print("\n--- Hypothetical Quantum-Gravitational Bridge ---")
    try:
        k_body_mass = float(input("Enter mass of the object in kg: "))
        max_distance_factor = float(input("Enter max distance as a multiple of Schwarzschild Radius (e.g., 5): "))
        schwarzschild_radius_val = calculate_schwarzschild_radius(k_body_mass)
        if schwarzschild_radius_val == 0:
            print("Invalid mass."); return
        distance_series = np.linspace(schwarzschild_radius_val + 1e-9, max_distance_factor * schwarzschild_radius_val, 200)
        explore_quantum_gravitational_differentiation(k_body_mass, distance_series, object_name="Massive Object")
    except ValueError:
        print("Invalid input. Please enter numerical values.")

def run_eos_stability_menu():
    print("\n--- EOS & Stability Explorer ---")
    sub = input("Choose: (N)ebula map or (C)osmic-strings? ").strip().lower()
    if sub == 'n':
        m_cloud_Msun = float(input("Cloud mass to compare (M_sun): "))
        mu = float(input("Mean molecular weight mu (e.g., 2.33): "))
        gamma = float(input("Adiabatic index gamma (e.g., 1.4 or 1.6667): "))
        rho_vals = np.logspace(-22, -16, 120)
        T_vals = np.logspace(1, 3.5, 120)
        plot_nebula_stability_grid(rho_vals, T_vals, mu=mu, gamma=gamma, mass_scale_Msun=m_cloud_Msun)
    elif sub == 'c':
        alpha_min = float(input("Wiggliness α min (e.g., 0): "))
        alpha_max = float(input("Wiggliness α max (e.g., 3): "))
        alphas = np.linspace(alpha_min, alpha_max, 200)
        plot_string_speeds_vs_wiggliness(alphas)
    else:
        print("Invalid choice.")

def run_database_query_menu():
    print("\n--- Astronomical Database Query ---")
    print("Select a database to query:")
    print("  1) ESA Gaia (by source ID)")
    print("  2) NEOSSAT (by coordinates)")
    print("  3) JAXA DARTS (by object name)")
    print("  4) NADC (by object name)")
    print("  5) STAC Index (by coordinates and time)")
    print("  6) ISRO (by object name)")
    try:
        db_choice = int(input("Enter your choice: "))
        if db_choice == 1:
            source_id = int(input("Enter a Gaia DR3 source ID (e.g., 5192257691883389568): "))
            fetch_gaia_data(source_id)
        elif db_choice == 2:
            ra_str = input("Enter Right Ascension (e.g., '10h40m10s'): ")
            dec_str = input("Enter Declination (e.g., '+22d30m5s'): ")
            coords = SkyCoord(ra=ra_str, dec=dec_str, frame='icrs')
            fetch_neossat_data(coords)
        elif db_choice == 3:
            target_name = input("Enter a celestial object name (e.g., 'M31'): ")
            fetch_darts_data(target_name)
        elif db_choice == 4:
            search_query = input("Enter a query (e.g., 'quasar'): ")
            fetch_nadc_data(search_query)
        elif db_choice == 5:
            ra = float(input("Enter RA in decimal degrees: "))
            dec = float(input("Enter Dec in decimal degrees: "))
            start_date = input("Enter start date (YYYY-MM-DD): ")
            end_date = input("Enter end date (YYYY-MM-DD): ")
            fetch_stac_data(ra, dec, start_date, end_date)
        elif db_choice == 6:
            target_name = input("Enter a celestial object name (e.g., 'AstroSat'): ")
            fetch_isro_data(target_name)
        else:
            print("Invalid choice.")
    except (ValueError, u.core.UnitConversionError) as e:
        print(f"Invalid input: {e}")

if __name__ == "__main__":
    print("--- Astrophysics & EOS Stability Explorer ---")
    print("This program explores several advanced physics concepts and databases.")
    while True:
        print("\n--- Main Menu ---")
        print("Select a topic to explore:")
        print("  1) Relativistic Effects on Matter & Radiation")
        print("  2) Hypothetical Quantum-Gravity Bridge")
        print("  3) Equation of State (EOS) & Stability")
        print("  4) International Astronomical Database Query")
        print("  5) Exit")
        main_choice = input("Enter your choice: ").strip().lower()

        if main_choice == '1':
            run_relativistic_comparison_menu()
        elif main_choice == '2':
            run_quantum_gravity_menu()
        elif main_choice == '3':
            run_eos_stability_menu()
        elif main_choice == '4':
            run_database_query_menu()
        elif main_choice == '5':
            print("\nExiting program. Thank you for exploring! 👋")
            break
        else:
            print("Invalid choice. Please enter a number from 1 to 5.")

Physics Simulator for finding EOS factors in Nebulae and Cosmic Strings

My love to ChatGPT and Google Gemini for the collaborative response efforts toward this astronomical exploration project!

import numpy as np
import matplotlib.pyplot as plt
import pubchempy as pcp
import requests
import json
import time
import re

# --- Scientific Constants ---
# Nuclear Physics Constants
CURIE_TO_BQ = 3.7e10           # Conversion factor: 1 Curie = 3.7 x 10^10 Becquerels (disintegrations per second)
LN2 = np.log(2)                # Natural logarithm of 2, used in half-life calculations

# Blackbody Radiation Constants
STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8 # Stefan-Boltzmann constant in W/(m^2·K^4)
PLANCK_CONSTANT = 6.62607015e-34           # Planck's constant in J·s
BOLTZMANN_CONSTANT = 1.380649e-23          # Boltzmann constant in J/K
SPEED_OF_LIGHT = 2.99792458e8              # Speed of light in a vacuum in m/s
WIEN_DISPLACEMENT_CONSTANT = 2.897771955e-3 # Wien's displacement law constant in m·K

# General Relativity Constants
GRAVITATIONAL_CONSTANT = 6.67430e-11 # G in m³·kg⁻¹·s⁻²
SPEED_OF_LIGHT_SQ = SPEED_OF_LIGHT**2

# --- Manual Nuclear Data for Gamma Emission ---
MANUAL_GAMMA_DATA = {
    "Cesium-137": {
        "half_life_days": 30.08 * 365.25,
        "gamma_energies_MeV": [0.662],
        "gamma_yields_per_decay": [0.851]
    },
    "Cobalt-60": {
        "half_life_days": 1925.21,
        "gamma_energies_MeV": [1.173, 1.332],
        "gamma_yields_per_decay": [0.999, 0.999]
    },
    "Iodine-131": {
        "half_life_days": 8.02,
        "gamma_energies_MeV": [0.364],
        "gamma_yields_per_decay": [0.817]
    },
    "Potassium-40": {
        "half_life_days": 1.251e9 * 365.25,
        "gamma_energies_MeV": [1.4608],
        "gamma_yields_per_decay": [0.1067]
    },
}

# ---------------------------------------------------------------------------
# PubChem Helper Functions
def _get_compound_cid_pubchem(compound_name):
    try:
        compounds = pcp.get_compounds(compound_name, 'name')
        return compounds[0].cid if compounds else None
    except Exception:
        return None

def _fetch_half_life_from_pubchem(cid):
    base_url = "https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/"
    url = f"{base_url}cid/{cid}/JSON/"
    try:
        response = requests.get(url, timeout=5)
        response.raise_for_status()
        data = response.json()
        search_sections = ['Atomic Mass, Half Life, and Decay', 'Chemical and Physical Properties', 'Depositor-Supplied']
        for record in data.get('Record', {}).get('Section', []):
            if record.get('TOCHeading') in search_sections:
                for prop_section in record.get('Section', []):
                    if 'half-life' in prop_section.get('TOCHeading', '').lower():
                        for item in prop_section.get('Information', []):
                            for value_dict in item.get('Value', {}).get('StringWithMarkup', []):
                                text = value_dict.get('String')
                                if text:
                                    match = re.match(r"(\d+\.?\d*(?:e[+-]?\d+)?)\s*(seconds|minutes|hours|days|years|s|min|hr|d|y)", text.lower())
                                    if match:
                                        return float(match.group(1)), match.group(2)
        return None, None
    except Exception:
        return None, None

def _convert_half_life_to_days(value, unit):
    if unit in ['s', 'seconds']: return value / (24 * 3600)
    if unit in ['min', 'minutes']: return value / (24 * 60)
    if unit in ['hr', 'hours']: return value / 24
    if unit in ['d', 'days']: return value
    if unit in ['y', 'years']: return value * 365.25
    return None

# ---------------------------------------------------------------------------
# Core Radiation Calculation Functions
def calculate_decay_constant(half_life_days):
    if half_life_days is None or half_life_days <= 0:
        return 0
    half_life_seconds = half_life_days * 24 * 3600
    return LN2 / half_life_seconds

def calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds):
    if decay_constant == 0:
        return initial_activity_Bq
    return initial_activity_Bq * np.exp(-decay_constant * time_elapsed_seconds)

def calculate_blackbody_spectral_radiance(wavelength_m, temperature_K):
    if temperature_K <= 0 or wavelength_m <= 0:
        return 0.0
    term1 = (2 * PLANCK_CONSTANT * SPEED_OF_LIGHT_SQ) / (wavelength_m**5)
    term2_exp_arg = (PLANCK_CONSTANT * SPEED_OF_LIGHT) / (wavelength_m * BOLTZMANN_CONSTANT * temperature_K)
    if term2_exp_arg > 700:
        return 0.0
    spectral_radiance = term1 / (np.exp(term2_exp_arg) - 1)
    return spectral_radiance

# ---------------------------------------------------------------------------
# General Relativity Functions
def calculate_schwarzschild_radius(mass_kg):
    """Calculates the Schwarzschild radius for a given mass."""
    if mass_kg <= 0:
        return 0.0
    return (2 * GRAVITATIONAL_CONSTANT * mass_kg) / SPEED_OF_LIGHT_SQ

def calculate_dilation_factor(mass_kg, distance_m):
    """
    Calculates the time dilation factor for an observer at infinity relative to an object in a gravitational field.
    A factor of 1 means no dilation. Values < 1 mean time is dilated (slowed down).
    """
    if mass_kg <= 0 or distance_m <= 0:
        return 1.0 # No dilation for zero mass or distance
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    if distance_m <= schwarzschild_radius:
        return 0.0 # At or inside event horizon, time stops for an outside observer
    return np.sqrt(1 - (schwarzschild_radius / distance_m))

def calculate_relativistic_activity(isotope_name, initial_activity_Bq, time_elapsed_days, mass_kg, distance_m):
    """
    Calculates activity as observed by a distant observer, accounting for gravitational time dilation
    and a hypothetical "time crystal" effect.
    """
    half_life_days = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("half_life_days")
    if half_life_days is None:
        cid = _get_compound_cid_pubchem(isotope_name)
        if cid:
            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)
            if hl_val:
                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)
    if half_life_days is None:
        return None

    time_elapsed_seconds = time_elapsed_days * 24 * 3600

    # Apply gravitational time dilation
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

    # An observer far away sees time running slower for the object, so the decay process appears slowed.
    effective_time_local_seconds = time_elapsed_seconds * dilation_factor

    # --- Hypothetical: curvature-based modifier ---
    curvature_parameter = calculate_gravitational_curvature_parameter(mass_kg, distance_m)
    time_crystal_modifier = _calculate_time_crystal_half_life_modifier(curvature_parameter)

    effective_half_life_days = half_life_days * time_crystal_modifier
    decay_constant = calculate_decay_constant(effective_half_life_days)
    current_activity_Bq = calculate_activity_at_time(initial_activity_Bq, decay_constant, effective_time_local_seconds)

    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("gamma_yields_per_decay", [0])
    total_gamma_yield_per_decay = sum(gamma_yields)
    total_gamma_disintegrations_per_second = current_activity_Bq * total_gamma_yield_per_decay
    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ
    return gamma_curies

def calculate_relativistic_peak_wavelength(local_temperature_K, mass_kg, distance_m):
    """Calculates the observed peak wavelength of a blackbody, accounting for gravitational redshift."""
    if local_temperature_K <= 0:
        return None
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    observed_temperature_K = local_temperature_K * dilation_factor
    if observed_temperature_K <= 0:
        return np.inf
    peak_wavelength_m = WIEN_DISPLACEMENT_CONSTANT / observed_temperature_K
    return peak_wavelength_m * 1e9 # Convert to nm

# ---------------------------------------------------------------------------
# Quantum and Gravitational Bridge Functions (hypothetical)
def calculate_gravitational_curvature_parameter(mass_kg, distance_m):
    """
    A simplified, unitless parameter representing local gravitational curvature.
    Higher values imply greater curvature.
    """
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    if distance_m <= 0 or distance_m <= schwarzschild_radius:
        return np.inf
    return schwarzschild_radius / distance_m

def quantum_differentiation_hypothetical(curvature_parameter, base_quantum_value):
    """
    Hypothetical function to illustrate 'quantum mechanical differentiation'.
    NOT based on established physics.
    """
    if curvature_parameter >= 1:
        return 0
    differentiated_value = base_quantum_value * (1 - curvature_parameter)
    return differentiated_value

def _calculate_time_crystal_half_life_modifier(curvature_parameter):
    """
    Hypothetical function to model a 'time crystal' effect driven by curvature.
    Returns a multiplicative modifier on the half-life.
    """
    if curvature_parameter >= 1:
        return 0.0
    modifier = np.tanh(1 - curvature_parameter)
    return modifier

# ---------------------------------------------------------------------------
# EOS & Stability (Nebulas + Cosmic Strings)
PROTON_MASS = 1.67262192369e-27
RADIATION_CONSTANT = 4 * STEFAN_BOLTZMANN_CONSTANT / SPEED_OF_LIGHT  # a_rad
K_B = BOLTZMANN_CONSTANT

# ---------- Gas / Nebula EOS ----------
def sound_speed_ideal_gas(T_K, mu=2.33, gamma=5/3):
    """
    Isentropic sound speed for an ideal gas:
      c_s = sqrt(gamma * k_B * T / (mu * m_p))
    mu ~ 2.33 for molecular clouds (H2 + He).
    """
    return np.sqrt(gamma * K_B * T_K / (mu * PROTON_MASS))

def pressure_ideal_gas(rho, T_K, mu=2.33):
    return rho * K_B * T_K / (mu * PROTON_MASS)

def pressure_polytrope(rho, K, gamma):
    return K * rho**gamma

def sound_speed_polytrope(rho, K, gamma):
    # c_s^2 = dP/dρ = gamma * K * ρ^(gamma-1)
    return np.sqrt(max(0.0, gamma * K * rho**(gamma - 1)))

def pressure_radiation(T_K):
    # P_rad = (1/3) a T^4
    return (1.0/3.0) * RADIATION_CONSTANT * T_K**4

def mixed_barotropic_pressure(rho, T_K, mu=2.33, gamma=5/3):
    # simple additive gas + radiation pressure
    return pressure_ideal_gas(rho, T_K, mu) + pressure_radiation(T_K)

# ---------- Nebula stability ----------
def jeans_length(rho, T_K, mu=2.33, gamma=5/3):
    """
    λ_J = c_s * sqrt(pi / (G * ρ))
    """
    cs = sound_speed_ideal_gas(T_K, mu, gamma)
    return cs * np.sqrt(np.pi / (GRAVITATIONAL_CONSTANT * rho))

def jeans_mass(rho, T_K, mu=2.33, gamma=5/3):
    """
    M_J ~ (4π/3) * (λ_J/2)^3 * ρ
    """
    lam = jeans_length(rho, T_K, mu, gamma)
    return (4*np.pi/3.0) * (lam/2.0)**3 * rho

def bonnor_ebert_mass(T_K, P_ext, mu=2.33):
    """
    Critical Bonnor–Ebert mass for isothermal sphere:
      M_BE ≈ 1.18 * (c_s^4) / (G^(3/2) * P_ext^(1/2))
    """
    cs = sound_speed_ideal_gas(T_K, mu, gamma=1.0)  # isothermal -> gamma=1
    return 1.18 * cs**4 / (GRAVITATIONAL_CONSTANT**1.5 * np.sqrt(P_ext))

# ---------- Cosmic strings (coarse-grained) ----------
def string_wiggliness_params(alpha=0.0, mu_base=1.0):
    """
    Phenomenology:
      U = μ (1+α),   T = μ /(1+α)
    α >= 0 encodes small-scale structure ("wiggles").
    """
    U = mu_base * (1.0 + alpha)
    T = mu_base / (1.0 + alpha)
    return U, T

def string_wave_speeds(alpha=0.0):
    """
    For the model above:
      c_T^2 = T/U = 1/(1+α)^2
    Use c_L ~ 1 (Nambu–Goto) as a simple default.
    """
    c_T = 1.0 / (1.0 + alpha)
    c_L = 1.0  # simple Nambu–Goto proxy
    return c_T, c_L

def string_effective_w(v_rms=0.5):
    """
    Effective EOS parameter for a network (toy model):
      w ≈ -1/3 * (1 - v_rms^2)
    """
    v2 = np.clip(v_rms**2, 0.0, 1.0)
    return -1.0/3.0 * (1.0 - v2)

def string_loop_oscillation_time(loop_length_m, alpha=0.0):
    """
    Fundamental oscillation time τ ~ L / (c_T * c).
    """
    c_T, _ = string_wave_speeds(alpha)
    return loop_length_m / (c_T * SPEED_OF_LIGHT)

# ---------- Visualizers ----------
def plot_nebula_stability_grid(rho_vals, T_vals, mu=2.33, gamma=5/3, mass_scale_Msun=1.0):
    """
    Heatmap: M_J / (mass_scale) to see where collapse is favored (M_J smaller than cloud mass).
    """
    from matplotlib.colors import LogNorm
    R, T = np.meshgrid(rho_vals, T_vals)
    MJ = np.vectorize(jeans_mass)(R, T, mu, gamma)
    plt.figure(figsize=(9,6))
    im = plt.pcolormesh(R, T, MJ/(mass_scale_Msun*1.98847e30), norm=LogNorm())
    plt.colorbar(im, label=r"$M_J/M_\odot$")
    plt.xscale('log'); plt.yscale('log')
    plt.xlabel(r"Density $\rho$ (kg m$^{-3}$)")
    plt.ylabel("Temperature (K)")
    plt.title("Nebula Stability Map (Jeans Mass in Solar Masses)")
    plt.tight_layout(); plt.show()

def plot_string_speeds_vs_wiggliness(alpha_vals):
    cT = [string_wave_speeds(a)[0] for a in alpha_vals]
    plt.figure(figsize=(8,5))
    plt.plot(alpha_vals, cT)
    plt.xlabel("Wiggliness α")
    plt.ylabel("Transverse wave speed c_T (units of c)")
    plt.title("Cosmic String Wave Speed vs Wiggliness")
    plt.grid(True); plt.tight_layout(); plt.show()

# ---------------------------------------------------------------------------
# Analysis Functions
def generate_normal_distribution_series(mean, std_dev, num_points=200, spread_factor=3):
    start = mean - spread_factor * std_dev
    end = mean + spread_factor * std_dev
    if start < 0:
        start = 0
    series_values = np.linspace(start, end, num_points)
    pdf_values = (1 / (std_dev * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((series_values - mean) / std_dev)**2)
    return series_values, pdf_values

def get_gamma_curies_single_point(isotope_name, initial_activity_Bq, time_elapsed_days):
    half_life_days = None
    if isotope_name in MANUAL_GAMMA_DATA:
        half_life_days = MANUAL_GAMMA_DATA[isotope_name]["half_life_days"]
    else:
        cid = _get_compound_cid_pubchem(isotope_name)
        if cid:
            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)
            if hl_val:
                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)
    if half_life_days is None:
        return None
    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("gamma_yields_per_decay", [0])
    total_gamma_yield_per_decay = sum(gamma_yields)
    decay_constant = calculate_decay_constant(half_life_days)
    time_elapsed_seconds = time_elapsed_days * 24 * 3600
    current_activity_Bq = calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds)
    total_gamma_disintegrations_per_second = current_activity_Bq * total_gamma_yield_per_decay
    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ
    return gamma_curies

def analyze_gamma_curies_comparison(isotope_name, initial_activity_Bq, mean_time_days, std_dev_time_days, mass_kg, distance_m, num_points=200):
    """
    Compares gamma curies decay in flat spacetime vs. near a massive k-body.
    """
    print(f"\n--- Comparing Gamma Curies: Flat Spacetime vs. Relativistic Environment ---")
    print(f"  Isotope: {isotope_name}, Initial Activity: {initial_activity_Bq:.2e} Bq")
    print(f"  Relativistic Environment: k-body Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m")
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")
    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")
    if dilation_factor < 0.99:
        print(f"  (A dilation factor less than 1 means time runs slower in the gravitational field.)")
    else:
        print(f"  (Dilation factor is near 1, indicating a weak gravitational effect.)")

    time_series_days, _ = generate_normal_distribution_series(mean_time_days, std_dev_time_days, num_points=num_points)

    # Flat Spacetime Calculation
    flat_curies_series = [get_gamma_curies_single_point(isotope_name, initial_activity_Bq, t_days) for t_days in time_series_days]

    # Relativistic Calculation
    relativistic_curies_series = [calculate_relativistic_activity(isotope_name, initial_activity_Bq, t_days, mass_kg, distance_m) for t_days in time_series_days]

    if all(c is None for c in flat_curies_series):
        print(f"  Calculation failed. Please ensure '{isotope_name}' has half-life data.")
        return

    # Plotting
    plt.figure(figsize=(12, 7))
    plt.plot(time_series_days, flat_curies_series, label='Flat Spacetime (Local Observer)')
    plt.plot(time_series_days, relativistic_curies_series, label=f'Relativistic (Distant Observer) w/ Curvature Modifier', linestyle='--')
    plt.title(f'Gamma Curies Decay: {isotope_name} - Flat vs. Relativistic')
    plt.xlabel('Time (days)')
    plt.ylabel('Gamma Curies (Ci)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

def simulate_blackbody_wien_comparison(temperatures_K, mass_kg, distance_m, wavelength_range_nm=(1, 10000)):
    """
    Compares blackbody peak wavelength in flat spacetime vs. near a massive k-body.
    """
    print(f"\n--- Comparing Blackbody Peak Wavelength: Flat vs. Relativistic ---")
    print(f"  Analyzing temperatures: {temperatures_K} K")
    print(f"  Relativistic Environment: k-body Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m")
    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)
    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)
    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")
    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")

    temp_list_sorted = sorted(temperatures_K)
    flat_peak_wl_list_nm = [(WIEN_DISPLACEMENT_CONSTANT / T) * 1e9 for T in temp_list_sorted if T > 0]
    relativistic_peak_wl_list_nm = [calculate_relativistic_peak_wavelength(T, mass_kg, distance_m) for T in temp_list_sorted if T > 0]

    # Plotting
    plt.figure(figsize=(8, 6))
    plt.plot(temp_list_sorted, flat_peak_wl_list_nm, 'o-', label='Flat Spacetime (Local Observer)')
    plt.plot(temp_list_sorted, relativistic_peak_wl_list_nm, 'r--', label='Relativistic (Distant Observer)')
    plt.title('Wien\'s Law: Peak Wavelength vs. Temperature (Flat vs. Relativistic)')
    plt.xlabel('Temperature (K)')
    plt.ylabel('Peak Wavelength (nm)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

def explore_quantum_gravitational_differentiation(mass_kg, distance_m_series):
    """
    New function to explore the hypothetical relationship between gravity and a quantum state.
    """
    print("\n--- Exploring a Hypothetical Gravitational-Quantum Bridge ---")
    print(f"  Analyzing Quantum State near a k-body with Mass={mass_kg:.2e} kg")
    curvature_series = [calculate_gravitational_curvature_parameter(mass_kg, d) for d in distance_m_series]
    base_quantum_value = 1.0
    differentiated_value_series = [quantum_differentiation_hypothetical(c, base_quantum_value) for c in curvature_series]

    plt.figure(figsize=(10, 6))
    plt.plot(distance_m_series, differentiated_value_series, 'b-')
    plt.axvline(x=calculate_schwarzschild_radius(mass_kg), color='r', linestyle='--', label='Schwarzschild Radius')
    plt.title('Hypothetical "Quantum Differentiation" vs. Distance from a k-body')
    plt.xlabel('Distance from k-body (m)')
    plt.ylabel('Hypothetical Differentiated Quantum Value (Unitless)')
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

# ---------------------------------------------------------------------------
# EOS Stability Explorer UI (Nebulas + Cosmic Strings)
def plot_nebula_stability_driver():
    try:
        m_cloud_Msun = float(input("  Cloud mass to compare (M_sun): "))
        mu = float(input("  Mean molecular weight mu (e.g., 2.33): "))
        gamma = float(input("  Adiabatic index gamma (e.g., 1.4 or 1.6667): "))
    except ValueError:
        print("  Invalid input.")
        return
    # Typical molecular cloud conditions in SI (kg/m^3) and K
    rho_vals = np.logspace(-22, -16, 120)  # kg/m^3
    T_vals = np.logspace(1, 3.5, 120)      # 10 K – ~3000 K
    plot_nebula_stability_grid(rho_vals, T_vals, mu=mu, gamma=gamma, mass_scale_Msun=m_cloud_Msun)
    print("  Darker regions (lower M_J/M_sun) imply easier collapse; where M_cloud > M_J, fragmentation likely.")

def plot_cosmic_string_driver():
    try:
        alpha_min = float(input("  Wiggliness α min (e.g., 0): "))
        alpha_max = float(input("  Wiggliness α max (e.g., 3): "))
    except ValueError:
        print("  Invalid input.")
        return
    alphas = np.linspace(alpha_min, alpha_max, 200)
    plot_string_speeds_vs_wiggliness(alphas)
    # quick demo of loop time
    try:
        L = float(input("  Example loop length L (meters): "))
        a_demo = float(input("  Example α (for τ): "))
    except ValueError:
        print("  Invalid input.")
        return
    tau = string_loop_oscillation_time(L, a_demo)
    print(f"  Loop oscillation time τ ≈ {tau:.3e} s at α={a_demo}. Lower c_T (higher α) → longer τ.")

# ---------------------------------------------------------------------------
# Main Program Execution
if __name__ == "__main__":
    print("--- Radiation + EOS Stability Explorer ---")
    print("\nThis program explores:")
    print("  1) Radioactive decay (flat vs relativistic, with a hypothetical curvature modifier).")
    print("  2) Blackbody redshift (flat vs relativistic).")
    print("  3) Hypothetical quantum–gravity differentiation.")
    print("  4) EOS Stability Explorer: Nebulas (Jeans/Bonnor–Ebert context) and Cosmic Strings (wiggliness & wave speeds).")

    while True:
        main_choice = input("Select: (R)elativistic Comparison, (Q)uantum-Gravity, (S)tability Explorer, (E)xit? ").strip().lower()
        if main_choice == 'e':
            break

        elif main_choice == 'r':
            sub_choice = input("Select Relativistic Scenario: (N)uclear Decay or (B)lackbody Radiation? ").strip().lower()
            try:
                # Get common relativistic parameters
                k_body_mass = float(input("\n  Enter mass of the k-body (e.g., a star) in kg: "))
                k_body_distance = float(input("  Enter distance from the center of the k-body in meters: "))
            except ValueError:
                print("  Invalid input. Please enter numerical values for mass and distance.")
                continue

            if sub_choice == 'n':
                print("\n--- Relativistic Nuclear Decay Comparison ---")
                isotope_name = input("  Enter radioactive isotope name (e.g., Cesium-137): ").strip()
                if isotope_name not in MANUAL_GAMMA_DATA and _get_compound_cid_pubchem(isotope_name) is None:
                    print(f"  '{isotope_name}' data not found.")
                    continue
                try:
                    initial_activity = float(input("  Enter initial activity in Bq: "))
                    mean_time = float(input("  Enter mean measurement time (days): "))
                    std_dev_time = float(input("  Enter standard deviation for time (days): "))
                except ValueError:
                    print("  Invalid input.")
                    continue
                analyze_gamma_curies_comparison(isotope_name, initial_activity, mean_time, std_dev_time, k_body_mass, k_body_distance)

            elif sub_choice == 'b':
                print("\n--- Relativistic Blackbody Radiation Comparison ---")
                temp_input_str = input("  Enter temperatures (K) separated by commas (e.g., 300, 1000, 5000): ")
                try:
                    temperatures_list = [float(t.strip()) for t in temp_input_str.split(',')]
                except ValueError:
                    print("  Invalid temperature input.")
                    continue
                simulate_blackbody_wien_comparison(temperatures_list, k_body_mass, k_body_distance)
            else:
                print("  Invalid choice. Please enter 'N' or 'B'.")

        elif main_choice == 'q':
            print("\n--- Quantum-Gravitational Bridge Exploration ---")
            try:
                k_body_mass = float(input("\n  Enter mass of the k-body (e.g., a black hole) in kg: "))
                max_distance_factor = float(input("  Enter max distance as a multiple of Schwarzschild Radius (e.g., 5): "))
                schwarzschild_radius_val = calculate_schwarzschild_radius(k_body_mass)
                if schwarzschild_radius_val == 0:
                    print("  Invalid mass. Cannot calculate Schwarzschild Radius.")
                    continue
                distance_series = np.linspace(schwarzschild_radius_val + 1e-9, max_distance_factor * schwarzschild_radius_val, 200)
                explore_quantum_gravitational_differentiation(k_body_mass, distance_series)
            except ValueError:
                print("  Invalid input. Please enter numerical values.")
                continue

        elif main_choice == 's':
            print("\n--- EOS Stability Explorer ---")
            sub = input("Choose: (N)ebula map, (C)osmic-strings? ").strip().lower()
            if sub == 'n':
                plot_nebula_stability_driver()
            elif sub == 'c':
                plot_cosmic_string_driver()
            else:
                print("  Invalid choice.")

        else:
            print("Invalid choice. Please enter 'R', 'Q', 'S', or 'E'.")

    print("\nExiting program. Thank you for exploring these concepts!")

Advanced Physics Simulator with Time Crystal Model

  Enjoy this Time-Crystal exploration program.

 

import numpy as np

import matplotlib.pyplot as plt

import pubchempy as pcp

import requests

import json

import time

import re

# --- Scientific Constants ---

# Nuclear Physics Constants

CURIE_TO_BQ = 3.7e10           # Conversion factor: 1 Curie = 3.7 x 10^10 Becquerels (disintegrations per second)

LN2 = np.log(2)                # Natural logarithm of 2, used in half-life calculations

# Blackbody Radiation Constants

STEFAN_BOLTZMANN_CONSTANT = 5.670374419e-8 # Stefan-Boltzmann constant in W/(m^2·K^4)

PLANCK_CONSTANT = 6.62607015e-34 # Planck's constant in J·s

BOLTZMANN_CONSTANT = 1.380649e-23 # Boltzmann constant in J/K

SPEED_OF_LIGHT = 2.99792458e8 # Speed of light in a vacuum in m/s

WIEN_DISPLACEMENT_CONSTANT = 2.897771955e-3 # Wien's displacement law constant in m·K

# General Relativity Constants (for the new functions)

GRAVITATIONAL_CONSTANT = 6.67430e-11 # G in m³·kg⁻¹·s⁻²

SPEED_OF_LIGHT_SQ = SPEED_OF_LIGHT**2

# --- Manual Nuclear Data for Gamma Emission ---

MANUAL_GAMMA_DATA = {

    "Cesium-137": {

        "half_life_days": 30.08 * 365.25,

        "gamma_energies_MeV": [0.662],

        "gamma_yields_per_decay": [0.851]

    },

    "Cobalt-60": {

        "half_life_days": 1925.21,

        "gamma_energies_MeV": [1.173, 1.332],

        "gamma_yields_per_decay": [0.999, 0.999]

    },

    "Iodine-131": {

        "half_life_days": 8.02,

        "gamma_energies_MeV": [0.364],

        "gamma_yields_per_decay": [0.817]

    },

    "Potassium-40": {

        "half_life_days": 1.251e9 * 365.25,

        "gamma_energies_MeV": [1.4608],

        "gamma_yields_per_decay": [0.1067]

    },

}

# --- PubChem Helper Functions (Unchanged) ---

def _get_compound_cid_pubchem(compound_name):

    try:

        compounds = pcp.get_compounds(compound_name, 'name')

        return compounds[0].cid if compounds else None

    except Exception:

        return None

def _fetch_half_life_from_pubchem(cid):

    base_url = "https://pubchem.ncbi.nlm.nih.gov/rest/pug_view/data/compound/"

    url = f"{base_url}cid/{cid}/JSON/"

    try:

        response = requests.get(url, timeout=5)

        response.raise_for_status()

        data = response.json()

        search_sections = ['Atomic Mass, Half Life, and Decay', 'Chemical and Physical Properties', 'Depositor-Supplied']

        for record in data.get('Record', {}).get('Section', []):

            if record.get('TOCHeading') in search_sections:

                for prop_section in record.get('Section', []):

                    if 'half-life' in prop_section.get('TOCHeading', '').lower():

                        for item in prop_section.get('Information', []):

                            for value_dict in item.get('Value', {}).get('StringWithMarkup', []):

                                text = value_dict.get('String')

                                if text:

                                    match = re.match(r"(\d+\.?\d*(?:e[+-]?\d+)?)\s*(seconds|minutes|hours|days|years|s|min|hr|d|y)", text.lower())

                                    if match:

                                        return float(match.group(1)), match.group(2)

        return None, None

    except Exception:

        return None, None

def _convert_half_life_to_days(value, unit):

    if unit in ['s', 'seconds']: return value / (24 * 3600)

    if unit in ['min', 'minutes']: return value / (24 * 60)

    if unit in ['hr', 'hours']: return value / 24

    if unit in ['d', 'days']: return value

    if unit in ['y', 'years']: return value * 365.25

    return None

# --- Core Radiation Calculation Functions (Modified to accept dilation factor) ---

def calculate_decay_constant(half_life_days):

    if half_life_days is None or half_life_days <= 0:

        return 0

    half_life_seconds = half_life_days * 24 * 3600

    return LN2 / half_life_seconds

def calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds):

    if decay_constant == 0:

        return initial_activity_Bq

    return initial_activity_Bq * np.exp(-decay_constant * time_elapsed_seconds)

def calculate_blackbody_spectral_radiance(wavelength_m, temperature_K):

    if temperature_K <= 0 or wavelength_m <= 0:

        return 0.0

    term1 = (2 * PLANCK_CONSTANT * SPEED_OF_LIGHT_SQ) / (wavelength_m**5)

    term2_exp_arg = (PLANCK_CONSTANT * SPEED_OF_LIGHT) / (wavelength_m * BOLTZMANN_CONSTANT * temperature_K)

    if term2_exp_arg > 700:

        return 0.0

    spectral_radiance = term1 / (np.exp(term2_exp_arg) - 1)

    return spectral_radiance

# --- UPDATED: General Relativity Functions ---

def calculate_schwarzschild_radius(mass_kg):

    """Calculates the Schwarzschild radius for a given mass."""

    if mass_kg <= 0:

        return 0.0

    return (2 * GRAVITATIONAL_CONSTANT * mass_kg) / SPEED_OF_LIGHT_SQ

def calculate_dilation_factor(mass_kg, distance_m):

    """

    Calculates the time dilation factor for an observer at infinity relative to an object in a gravitational field.

    A factor of 1 means no dilation. Values < 1 mean time is dilated (slowed down).

    """

    if mass_kg <= 0 or distance_m <= 0:

        return 1.0 # No dilation for zero mass or distance

   

    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)

    if distance_m <= schwarzschild_radius:

        return 0.0 # At or inside event horizon, time stops for an outside observer

   

    return np.sqrt(1 - (schwarzschild_radius / distance_m))

def calculate_relativistic_activity(isotope_name, initial_activity_Bq, time_elapsed_days, mass_kg, distance_m):

    """

    Calculates activity as observed by a distant observer, accounting for gravitational time dilation

    and a hypothetical "time crystal" effect.

    """

    half_life_days = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("half_life_days")

    if half_life_days is None:

        cid = _get_compound_cid_pubchem(isotope_name)

        if cid:

            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)

            if hl_val:

                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)

    if half_life_days is None:

        return None

    time_elapsed_seconds = time_elapsed_days * 24 * 3600

    # Apply gravitational time dilation

    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

   

    # An observer far away sees time running slower for the object, so the decay process appears slowed.

    effective_time_local_seconds = time_elapsed_seconds * dilation_factor

    # --- NEW: Apply Time Crystal "dynamic range modifier" ---

    # We use a hypothetical model where the half-life itself is affected by extreme gravity.

    # The `curvature_parameter` serves as a proxy for the 'valence instability' and the

    # potential for the system to "collimate" into a time crystal state.

    curvature_parameter = calculate_gravitational_curvature_parameter(mass_kg, distance_m)

    time_crystal_modifier = _calculate_time_crystal_half_life_modifier(curvature_parameter)

   

    # We multiply the original half-life by this modifier. A modifier < 1 means the

    # "effective" half-life shortens, speeding up the perceived decay rate, which

    # could model the collapse into a new, more stable state.

    effective_half_life_days = half_life_days * time_crystal_modifier

   

    decay_constant = calculate_decay_constant(effective_half_life_days)

   

    current_activity_Bq = calculate_activity_at_time(initial_activity_Bq, decay_constant, effective_time_local_seconds)

   

    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("gamma_yields_per_decay", [0])

    total_gamma_yield_per_decay = sum(gamma_yields)

    total_gamma_disintegrations_per_second = current_activity_Bq * total_gamma_yield_per_decay

    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ

    return gamma_curies

def calculate_relativistic_peak_wavelength(local_temperature_K, mass_kg, distance_m):

    """

    Calculates the observed peak wavelength of a blackbody, accounting for gravitational redshift.

    """

    if local_temperature_K <= 0:

        return None

    # Gravitational redshift means a distant observer sees the blackbody as cooler.

    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

    observed_temperature_K = local_temperature_K * dilation_factor

    # Use Wien's law with the observed temperature

    if observed_temperature_K <= 0:

        return np.inf # Effectively no peak wavelength observed

   

    peak_wavelength_m = WIEN_DISPLACEMENT_CONSTANT / observed_temperature_K

    return peak_wavelength_m * 1e9 # Convert to nm

# --- NEW: Quantum and Gravitational Bridge Functions ---

def calculate_gravitational_curvature_parameter(mass_kg, distance_m):

    """

    A simplified, unitless parameter representing local gravitational curvature.

    Higher values imply greater curvature.

    """

    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)

    if distance_m <= 0 or distance_m <= schwarzschild_radius:

        return np.inf

    return schwarzschild_radius / distance_m

def quantum_differentiation_hypothetical(curvature_parameter, base_quantum_value):

    """

    Hypothetical function to illustrate 'quantum mechanical differentiation'.

    This is NOT based on established physics but serves as an illustrative model.

    It proposes that a quantum value (e.g., a spin state) changes with gravitational curvature.

   

    For this example, we'll model it as a simple linear relationship for visualization.

    """

    if curvature_parameter >= 1:

        return 0 # At or inside the event horizon, the effect is undefined or at its maximum

   

    # A hypothetical effect: the base value is 'differentiated' by the curvature.

    # We use (1 - curvature_parameter) to get a value that goes from 1 to 0 as curvature increases.

    differentiated_value = base_quantum_value * (1 - curvature_parameter)

    return differentiated_value

def _calculate_time_crystal_half_life_modifier(curvature_parameter):

    """

    Hypothetical function to model the "time crystal" effect.

    As gravitational curvature increases, the half-life is modified, simulating a

    "collapsing" or "colliding" into a new, stable-like state.

   

    This is not based on established physics.

    The function returns a modifier that, when multiplied by the half-life,

    changes its value. A value < 1 shortens the half-life.

    """

    if curvature_parameter >= 1:

        # At or inside the event horizon, the half-life effectively becomes zero.

        # This models the ultimate collapse into a "time crystal" state.

        return 0.0

    # We use a simple, non-linear function to model the effect.

    # As curvature approaches 1, the modifier approaches zero.

    modifier = np.tanh(1 - curvature_parameter)

    return modifier

# --- Utility Functions (Unchanged) ---

def generate_normal_distribution_series(mean, std_dev, num_points=200, spread_factor=3):

    start = mean - spread_factor * std_dev

    end = mean + spread_factor * std_dev

    if start < 0:

        start = 0

    series_values = np.linspace(start, end, num_points)

    pdf_values = (1 / (std_dev * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((series_values - mean) / std_dev)**2)

    return series_values, pdf_values

def get_gamma_curies_single_point(isotope_name, initial_activity_Bq, time_elapsed_days):

    half_life_days = None

    if isotope_name in MANUAL_GAMMA_DATA:

        half_life_days = MANUAL_GAMMA_DATA[isotope_name]["half_life_days"]

    else:

        cid = _get_compound_cid_pubchem(isotope_name)

        if cid:

            hl_val, hl_unit = _fetch_half_life_from_pubchem(cid)

            if hl_val:

                half_life_days = _convert_half_life_to_days(hl_val, hl_unit)

    if half_life_days is None:

        return None

    gamma_yields = MANUAL_GAMMA_DATA.get(isotope_name, {}).get("gamma_yields_per_decay", [0])

    total_gamma_yield_per_decay = sum(gamma_yields)

    decay_constant = calculate_decay_constant(half_life_days)

    time_elapsed_seconds = time_elapsed_days * 24 * 3600

    current_activity_Bq = calculate_activity_at_time(initial_activity_Bq, decay_constant, time_elapsed_seconds)

    total_gamma_disintegrations_per_second = current_activity_Bq * total_gamma_yield_per_decay

    gamma_curies = total_gamma_disintegrations_per_second / CURIE_TO_BQ

    return gamma_curies

# --- Analysis Functions (Modified to include new relativistic options) ---

def analyze_gamma_curies_comparison(isotope_name, initial_activity_Bq, mean_time_days, std_dev_time_days, mass_kg, distance_m, num_points=200):

    """

    Compares gamma curies decay in flat spacetime vs. near a massive k-body.

    """

    print(f"\n--- Comparing Gamma Curies: Flat Spacetime vs. Relativistic Environment ---")

    print(f"  Isotope: {isotope_name}, Initial Activity: {initial_activity_Bq:.2e} Bq")

    print(f"  Relativistic Environment: k-body Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m")

   

    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)

    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")

    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")

    if dilation_factor < 0.99:

        print(f"  (A dilation factor less than 1 means time runs slower in the gravitational field.)")

    else:

        print(f"  (Dilation factor is near 1, indicating a weak gravitational effect.)")

    time_series_days, _ = generate_normal_distribution_series(mean_time_days, std_dev_time_days, num_points=num_points)

    # Flat Spacetime Calculation

    flat_curies_series = [get_gamma_curies_single_point(isotope_name, initial_activity_Bq, t_days) for t_days in time_series_days]

    # Relativistic Calculation

    relativistic_curies_series = [calculate_relativistic_activity(isotope_name, initial_activity_Bq, t_days, mass_kg, distance_m) for t_days in time_series_days]

    if all(c is None for c in flat_curies_series):

        print(f"  Calculation failed. Please ensure '{isotope_name}' has half-life data.")

        return

    # Plotting

    plt.figure(figsize=(12, 7))

    plt.plot(time_series_days, flat_curies_series, label='Flat Spacetime (Local Observer)')

    plt.plot(time_series_days, relativistic_curies_series, label=f'Relativistic (Distant Observer) w/ Time Crystal Modifier', linestyle='--')

    plt.title(f'Gamma Curies Decay: {isotope_name} - Flat vs. Relativistic')

    plt.xlabel('Time (days)')

    plt.ylabel('Gamma Curies (Ci)')

    plt.grid(True)

    plt.legend()

    plt.tight_layout()

    plt.show()

def simulate_blackbody_wien_comparison(temperatures_K, mass_kg, distance_m, wavelength_range_nm=(1, 10000)):

    """

    Compares blackbody peak wavelength in flat spacetime vs. near a massive k-body.

    """

    print(f"\n--- Comparing Blackbody Peak Wavelength: Flat vs. Relativistic ---")

    print(f"  Analyzing temperatures: {temperatures_K} K")

    print(f"  Relativistic Environment: k-body Mass={mass_kg:.2e} kg, Distance={distance_m:.2e} m")

    schwarzschild_radius = calculate_schwarzschild_radius(mass_kg)

    dilation_factor = calculate_dilation_factor(mass_kg, distance_m)

    print(f"  Schwarzschild Radius: {schwarzschild_radius:.2e} m")

    print(f"  Gravitational Time Dilation Factor: {dilation_factor:.4f}")

    temp_list_sorted = sorted(temperatures_K)

    flat_peak_wl_list_nm = [(WIEN_DISPLACEMENT_CONSTANT / T) * 1e9 for T in temp_list_sorted if T > 0]

   

    relativistic_peak_wl_list_nm = [calculate_relativistic_peak_wavelength(T, mass_kg, distance_m) for T in temp_list_sorted if T > 0]

    # Plotting

    plt.figure(figsize=(8, 6))

    plt.plot(temp_list_sorted, flat_peak_wl_list_nm, 'o-', label='Flat Spacetime (Local Observer)')

    plt.plot(temp_list_sorted, relativistic_peak_wl_list_nm, 'r--', label='Relativistic (Distant Observer)')

    plt.title('Wien\'s Law: Peak Wavelength vs. Temperature (Flat vs. Relativistic)')

    plt.xlabel('Temperature (K)')

    plt.ylabel('Peak Wavelength (nm)')

    plt.grid(True)

    plt.legend()

    plt.tight_layout()

    plt.show()

def explore_quantum_gravitational_differentiation(mass_kg, distance_m_series):

    """

    New function to explore the hypothetical relationship between gravity and a quantum state.

    """

    print("\n--- Exploring a Hypothetical Gravitational-Quantum Bridge ---")

    print(f"  Analyzing Quantum State near a k-body with Mass={mass_kg:.2e} kg")

    curvature_series = [calculate_gravitational_curvature_parameter(mass_kg, d) for d in distance_m_series]

   

    # Let's use a base quantum value, for example, a spin state (up/down).

    # We'll use 1.0 as a base value to represent an "undifferentiated" quantum state.

    base_quantum_value = 1.0

   

    differentiated_value_series = [quantum_differentiation_hypothetical(c, base_quantum_value) for c in curvature_series]

    plt.figure(figsize=(10, 6))

    plt.plot(distance_m_series, differentiated_value_series, 'b-')

    plt.axvline(x=calculate_schwarzschild_radius(mass_kg), color='r', linestyle='--', label='Schwarzschild Radius')

    plt.title('Hypothetical "Quantum Differentiation" vs. Distance from a k-body')

    plt.xlabel('Distance from k-body (m)')

    plt.ylabel('Hypothetical Differentiated Quantum Value (Unitless)')

    plt.grid(True)

    plt.legend()

    plt.tight_layout()

    plt.show()

# --- Main Program Execution (Modified) ---

if __name__ == "__main__":

    print("--- Radiation Analysis Program: Interpreting Data Behavior ---")

    print("\nThis program explores two physics concepts in both a standard, flat-spacetime environment and a relativistic one.")

   

    print("1. **Radioactive Decay:** Compares the decay rate of an isotope as seen by a local observer vs. a distant one in a gravitational field, including a hypothetical 'time crystal' effect.")

    print("2. **Blackbody Radiation:** Compares the peak emission wavelength of a blackbody as seen by a local observer vs. a distant one (gravitational redshift).")

    print("3. **Hypothetical Quantum-Gravity Bridge:** A new, conceptual section to explore how a quantum state might be 'differentiated' by a gravitational field.\n")

    while True:

        main_choice = input("Select: (R)elativistic Comparison, (Q)uantum-Gravity, (E)xit? ").strip().lower()

        if main_choice == 'e':

            break

        elif main_choice == 'r':

            sub_choice = input("Select Relativistic Scenario: (N)uclear Decay or (B)lackbody Radiation? ").strip().lower()

            try:

                # Get common relativistic parameters

                k_body_mass = float(input("\n  Enter mass of the k-body (e.g., a star) in kg: "))

                k_body_distance = float(input("  Enter distance from the center of the k-body in meters: "))

            except ValueError:

                print("  Invalid input. Please enter numerical values for mass and distance.")

                continue

            if sub_choice == 'n':

                print("\n--- Relativistic Nuclear Decay Comparison ---")

                isotope_name = input("  Enter radioactive isotope name (e.g., Cesium-137): ").strip()

                if isotope_name not in MANUAL_GAMMA_DATA and _get_compound_cid_pubchem(isotope_name) is None:

                    print(f"  '{isotope_name}' data not found.")

                    continue

                try:

                    initial_activity = float(input("  Enter initial activity in Bq: "))

                    mean_time = float(input("  Enter mean measurement time (days): "))

                    std_dev_time = float(input("  Enter standard deviation for time (days): "))

                except ValueError:

                    print("  Invalid input.")

                    continue

               

                analyze_gamma_curies_comparison(isotope_name, initial_activity, mean_time, std_dev_time, k_body_mass, k_body_distance)

            elif sub_choice == 'b':

                print("\n--- Relativistic Blackbody Radiation Comparison ---")

                temp_input_str = input("  Enter temperatures (K) separated by commas (e.g., 300, 1000, 5000): ")

                try:

                    temperatures_list = [float(t.strip()) for t in temp_input_str.split(',')]

                except ValueError:

                    print("  Invalid temperature input.")

                    continue

                simulate_blackbody_wien_comparison(temperatures_list, k_body_mass, k_body_distance)

            else:

                print("  Invalid choice. Please enter 'N' or 'B'.")

        elif main_choice == 'q':

            print("\n--- Quantum-Gravitational Bridge Exploration ---")

            try:

                k_body_mass = float(input("\n  Enter mass of the k-body (e.g., a black hole) in kg: "))

                # We'll generate a range of distances to show the effect

                max_distance_factor = float(input("  Enter max distance as a multiple of Schwarzschild Radius (e.g., 5): "))

               

                schwarzschild_radius_val = calculate_schwarzschild_radius(k_body_mass)

                if schwarzschild_radius_val == 0:

                    print("  Invalid mass. Cannot calculate Schwarzschild Radius.")

                    continue

               

                distance_series = np.linspace(schwarzschild_radius_val + 1e-9, max_distance_factor * schwarzschild_radius_val, 200)

                explore_quantum_gravitational_differentiation(k_body_mass, distance_series)

            except ValueError:

                print("  Invalid input. Please enter numerical values.")

                continue

        else:

            print("Invalid choice. Please enter 'R', 'Q', or 'E'.")

    print("\nExiting program. Thank you for exploring these concepts!")