Update PySDM/mcp_output/mcp_plugin/mcp_service.py

PySDM/mcp_output/mcp_plugin/mcp_service.py

@@ -1,87 +1,400 @@
+import os
+import sys
+from typing import List, Optional, Dict, Any
+import math
+
+# Add the local source directory to sys.path
+source_path = os.path.join(os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))), "source")
+if source_path not in sys.path:
+    sys.path.insert(0, source_path)
+
 from fastmcp import FastMCP
+import numpy as np
+
+# Import PySDM modules
+from PySDM.physics import constants as const
+from PySDM.physics import si
+from PySDM.physics.trivia import Trivia
+from PySDM.formulae import Formulae
 
 # Create the FastMCP service application
 mcp = FastMCP("pysdm_service")
 
 
-@mcp.tool(name="list_attributes", description="List all available attributes in PySDM")
-def list_attributes() -> dict:
+# ===================== Physical Constants =====================
+
+@mcp.tool(name="get_physical_constants", description="Retrieve physical constants")
+def get_physical_constants() -> dict:
+    """
+    Retrieve physical constants used in PySDM simulations.
+
+    Returns:
+    - dict: Dictionary containing physical constants with their values and units.
+    """
+    try:
+        result = {
+            "fundamental_constants": {
+                "R_str": {"value": float(const.sci.R), "unit": "J/(K·mol)", "description": "Universal gas constant"},
+                "N_A": {"value": float(const.sci.N_A), "unit": "1/mol", "description": "Avogadro constant"},
+                "g_std": {"value": float(const.sci.g), "unit": "m/s²", "description": "Standard gravity"},
+                "PI": {"value": float(const.PI), "unit": "dimensionless", "description": "Pi"},
+            },
+            "thermodynamic_constants": {
+                "T0": {"value": float(const.T0 / si.kelvin), "unit": "K", "description": "Zero Celsius in Kelvin"},
+                "sqrt_two": {"value": float(const.sqrt_two), "unit": "dimensionless", "description": "Square root of 2"},
+                "sqrt_pi": {"value": float(const.sqrt_pi), "unit": "dimensionless", "description": "Square root of pi"},
+            },
+            "numerical_constants": {
+                "ONE_THIRD": {"value": float(const.ONE_THIRD), "unit": "dimensionless", "description": "1/3"},
+                "TWO_THIRDS": {"value": float(const.TWO_THIRDS), "unit": "dimensionless", "description": "2/3"},
+                "PI_4_3": {"value": float(const.PI_4_3), "unit": "dimensionless", "description": "4π/3"},
+            },
+            "concentration_units": {
+                "PPM": {"value": float(const.PPM), "unit": "dimensionless", "description": "Parts per million"},
+                "PPB": {"value": float(const.PPB), "unit": "dimensionless", "description": "Parts per billion"},
+                "PER_CENT": {"value": float(const.PER_CENT), "unit": "dimensionless", "description": "Percent"},
+                "PER_MILLE": {"value": float(const.PER_MILLE), "unit": "dimensionless", "description": "Per mille"},
+            }
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Formulae Tools =====================
+
+@mcp.tool(name="calculate_saturation_vapour_pressure", description="Calculate saturation vapour pressure over water")
+def calculate_saturation_vapour_pressure(temperature_kelvin: float, method: str = "FlatauWalkoCotton") -> dict:
+    """
+    Calculate saturation vapour pressure over liquid water using PySDM Formulae.
+
+    Parameters:
+    - temperature_kelvin (float): Temperature in Kelvin.
+    - method (str): Method to use - 'FlatauWalkoCotton', 'AugustRocheMagnus', 'MurphyKoop2005', etc.
+
+    Returns:
+    - dict: Saturation vapour pressure in Pascals and hectopascals.
+    """
+    try:
+        formulae = Formulae(saturation_vapour_pressure=method)
+        T = temperature_kelvin * si.kelvin
+        pvs = formulae.saturation_vapour_pressure.pvs_water(formulae.constants, T)
+        pvs_value = float(pvs / si.pascal)
+        
+        result = {
+            "temperature_K": temperature_kelvin,
+            "temperature_C": temperature_kelvin - 273.15,
+            "saturation_vapour_pressure_Pa": pvs_value,
+            "saturation_vapour_pressure_hPa": pvs_value / 100,
+            "method": method
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+@mcp.tool(name="calculate_condensation", description="Calculate condensation rates")
+def calculate_condensation(temperature: float, pressure: float, relative_humidity: float = 1.0) -> dict:
+    """
+    Calculate condensation-related parameters using PySDM.
+
+    Parameters:
+    - temperature (float): Temperature in Kelvin.
+    - pressure (float): Pressure in Pascals.
+    - relative_humidity (float): Relative humidity (0-1 or as fraction >1 for supersaturation).
+
+    Returns:
+    - dict: Condensation parameters including supersaturation and vapour pressure.
+    """
+    try:
+        formulae = Formulae()
+        T = temperature * si.kelvin
+        pvs = formulae.saturation_vapour_pressure.pvs_water(formulae.constants, T)
+        pvs_value = float(pvs / si.pascal)
+        
+        # Actual vapour pressure
+        pv = relative_humidity * pvs_value
+        
+        # Supersaturation
+        supersaturation = relative_humidity - 1.0
+        
+        # Water vapour mixing ratio (eps = Mv/Md ≈ 0.622)
+        eps = 0.622
+        mixing_ratio = eps * pv / (pressure - pv) if pressure > pv else float('nan')
+        
+        # Specific humidity
+        specific_humidity = mixing_ratio / (1 + mixing_ratio) if not np.isnan(mixing_ratio) else float('nan')
+        
+        result = {
+            "temperature_K": temperature,
+            "pressure_Pa": pressure,
+            "saturation_vapour_pressure_Pa": pvs_value,
+            "actual_vapour_pressure_Pa": pv,
+            "relative_humidity": relative_humidity,
+            "supersaturation": supersaturation,
+            "supersaturation_percent": supersaturation * 100,
+            "water_vapour_mixing_ratio": mixing_ratio,
+            "specific_humidity": specific_humidity
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Particle Dynamics =====================
+
+@mcp.tool(name="simulate_particles", description="Simulate particle dynamics using PySDM")
+def simulate_particles(particle_count: int, time_step: float) -> dict:
     """
-    List all available attributes in PySDM.
+    Get information about particle simulation parameters in PySDM.
+
+    Parameters:
+    - particle_count (int): Number of super-droplets to simulate.
+    - time_step (float): Time step for the simulation in seconds.
 
     Returns:
-    - dict: A dictionary with success status and list of available attributes.
+    - dict: Simulation configuration and recommendations.
     """
     try:
-        from PySDM.attributes import __all__ as attributes
-        return {
-            "success": True,
-            "attributes": attributes,
-            "count": len(attributes)
+        from PySDM.dynamics.condensation import DEFAULTS as COND_DEFAULTS
+        
+        defaults = {
+            "rtol_x": COND_DEFAULTS.rtol_x,
+            "rtol_thd": COND_DEFAULTS.rtol_thd,
+            "dt_cond_range": (float(COND_DEFAULTS.cond_range[0] / si.second), float(COND_DEFAULTS.cond_range[1] / si.second)),
+            "schedule": COND_DEFAULTS.schedule,
+        }
+        
+        result = {
+            "configuration": {
+                "n_sd": particle_count,
+                "dt": time_step,
+                "dt_unit": "seconds"
+            },
+            "solver_defaults": defaults,
+            "available_dynamics": [
+                "Condensation",
+                "Collision/Coalescence",
+                "Displacement",
+                "Freezing",
+                "AqueousChemistry",
+                "IsotopicFractionation",
+                "VapourDepositionOnIce"
+            ],
+            "recommendations": {
+                "adaptive_timestep": "Recommended for condensation",
+                "suggested_n_sd": "100-10000 for typical cloud simulations"
+            }
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
 
-@mcp.tool(name="list_dynamics", description="List all available dynamics in PySDM")
-def list_dynamics() -> dict:
+# ===================== Trivia/Utility Functions =====================
+
+@mcp.tool(name="calculate_droplet_volume", description="Calculate droplet volume from radius")
+def calculate_droplet_volume(radius_um: float) -> dict:
     """
-    List all available dynamics in PySDM.
+    Calculate droplet volume and related properties using PySDM Trivia functions.
+
+    Parameters:
+    - radius_um (float): Droplet radius in micrometers.
 
     Returns:
-    - dict: A dictionary with success status and list of available dynamics.
+    - dict: Volume, surface area, and mass of the droplet.
     """
     try:
-        from PySDM.dynamics import __all__ as dynamics
-        return {
-            "success": True,
-            "dynamics": dynamics,
-            "count": len(dynamics)
+        formulae = Formulae()
+        radius_m = radius_um * 1e-6
+        
+        volume = Trivia.volume(formulae.constants, radius_m)
+        surface_area = Trivia.area(formulae.constants, radius_m)
+        
+        # Assuming water density ~1000 kg/m³
+        rho_w = 1000.0
+        mass = rho_w * float(volume)
+        
+        result = {
+            "radius_um": radius_um,
+            "radius_m": radius_m,
+            "volume_m3": float(volume),
+            "volume_um3": float(volume) * 1e18,
+            "surface_area_m2": float(surface_area),
+            "surface_area_um2": float(surface_area) * 1e12,
+            "mass_kg": mass,
+            "mass_ng": mass * 1e12
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
 
-@mcp.tool(name="list_physics", description="List all available physics modules in PySDM")
-def list_physics() -> dict:
+@mcp.tool(name="calculate_radius_from_volume", description="Calculate droplet radius from volume")
+def calculate_radius_from_volume(volume_um3: float) -> dict:
     """
-    List all available physics modules in PySDM.
+    Calculate droplet radius from volume using PySDM Trivia functions.
+
+    Parameters:
+    - volume_um3 (float): Droplet volume in cubic micrometers.
 
     Returns:
-    - dict: A dictionary with success status and list of available physics modules.
+    - dict: Radius in various units.
     """
     try:
-        from PySDM.physics import __all__ as physics_modules
-        return {
-            "success": True,
-            "physics_modules": physics_modules,
-            "count": len(physics_modules)
+        formulae = Formulae()
+        volume_m3 = volume_um3 * 1e-18
+        
+        radius_m = Trivia.radius(formulae.constants, volume_m3)
+        radius_um = float(radius_m) * 1e6
+        
+        result = {
+            "volume_um3": volume_um3,
+            "volume_m3": volume_m3,
+            "radius_m": float(radius_m),
+            "radius_um": radius_um,
+            "diameter_um": 2 * radius_um
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
+
 
+# ===================== Kappa-Köhler Hygroscopicity =====================
 
-@mcp.tool(name="simulate", description="Run a PySDM simulation")
-def simulate(config: dict) -> dict:
+@mcp.tool(name="calculate_kappa_koehler", description="Calculate critical supersaturation using kappa-Köhler theory")
+def calculate_kappa_koehler(dry_radius_um: float, kappa: float, temperature_kelvin: float = 293.15) -> dict:
     """
-    Run a PySDM simulation based on the provided configuration.
+    Calculate critical supersaturation and radius using kappa-Köhler theory.
 
     Parameters:
-    - config: A dictionary containing simulation parameters.
+    - dry_radius_um (float): Dry aerosol radius in micrometers.
+    - kappa (float): Hygroscopicity parameter (kappa).
+    - temperature_kelvin (float): Temperature in Kelvin (default 293.15 K = 20°C).
+
+    Returns:
+    - dict: Critical supersaturation and activation radius.
+    """
+    try:
+        formulae = Formulae(hygroscopicity="KappaKoehlerLeadingTerms")
+        
+        # Physical constants
+        sigma = 0.072  # Surface tension of water (N/m)
+        Mv = 18.015e-3  # kg/mol
+        rho_w = 1000.0  # kg/m³
+        R = const.sci.R
+        T = temperature_kelvin
+        
+        dry_radius_m = dry_radius_um * 1e-6
+        
+        # Kelvin parameter A
+        A = 2 * sigma * Mv / (rho_w * R * T)
+        
+        # Critical supersaturation (approximation from leading terms)
+        S_c = math.sqrt(4 * A**3 / (27 * kappa * dry_radius_m**3))
+        
+        # Critical radius
+        r_c = math.sqrt(3 * kappa * dry_radius_m**3 / A)
+        
+        result = {
+            "dry_radius_um": dry_radius_um,
+            "kappa": kappa,
+            "temperature_K": temperature_kelvin,
+            "kelvin_parameter_A": A,
+            "critical_supersaturation": S_c,
+            "critical_supersaturation_percent": S_c * 100,
+            "critical_radius_um": r_c * 1e6,
+            "activation_diameter_um": 2 * r_c * 1e6
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Isotope Tools =====================
+
+@mcp.tool(name="get_isotope_constants", description="Get water isotope constants")
+def get_isotope_constants() -> dict:
+    """
+    Get water isotope constants (VSMOW standard) from PySDM.
+
+    Returns:
+    - dict: Isotope abundance ratios and atomic masses.
+    """
+    try:
+        from PySDM.physics import constants_defaults as cd
+        
+        result = {
+            "VSMOW_ratios": {
+                "R_2H": {"value": float(cd.VSMOW_R_2H), "description": "Deuterium abundance ratio"},
+                "R_3H": {"value": float(cd.VSMOW_R_3H), "description": "Tritium abundance ratio"},
+                "R_18O": {"value": float(cd.VSMOW_R_18O), "description": "Oxygen-18 abundance ratio"},
+                "R_17O": {"value": float(cd.VSMOW_R_17O), "description": "Oxygen-17 abundance ratio"},
+            },
+            "atomic_masses_kg_per_mol": {
+                "M_1H": float(cd.M_1H / (si.g / si.mole)),
+                "M_2H": float(cd.M_2H / (si.g / si.mole)),
+                "M_16O": float(cd.M_16O / (si.g / si.mole)),
+                "M_18O": float(cd.M_18O / (si.g / si.mole)),
+            },
+            "description": "VSMOW (Vienna Standard Mean Ocean Water) is the international standard for water isotope ratios"
+        }
+        return {"success": True, "result": result, "error": None}
+    except Exception as e:
+        return {"success": False, "result": None, "error": str(e)}
+
+
+# ===================== Available Formulae =====================
+
+@mcp.tool(name="list_available_formulae", description="List available physics formulae in PySDM")
+def list_available_formulae() -> dict:
+    """
+    List available physics formulae options in PySDM.
 
     Returns:
-    - dict: Simulation results or status.
+    - dict: Categories of formulae with available options.
     """
     try:
-        from PySDM import Simulation
-        simulation = Simulation(config)
-        results = simulation.run()
-        return {
-            "success": True,
-            "results": results
+        result = {
+            "saturation_vapour_pressure": [
+                "FlatauWalkoCotton",
+                "AugustRocheMagnus",
+                "Lowe1977",
+                "MurphyKoop2005",
+                "Wexler1976",
+                "Bolton1980"
+            ],
+            "hygroscopicity": [
+                "KappaKoehler",
+                "KappaKoehlerLeadingTerms"
+            ],
+            "latent_heat_vapourisation": [
+                "Kirchhoff",
+                "Constant"
+            ],
+            "drop_growth": [
+                "Mason1971",
+                "FuchsSutugin"
+            ],
+            "surface_tension": [
+                "Constant",
+                "CompressedFilm"
+            ],
+            "terminal_velocity": [
+                "GunnKinzer1949",
+                "PowerSeries",
+                "RogersYau"
+            ],
+            "freezing_temperature_spectrum": [
+                "Null",
+                "Bigg1953",
+                "Niemand_et_al_2012"
+            ],
+            "description": "These are configurable physics options in PySDM.Formulae"
         }
+        return {"success": True, "result": result, "error": None}
     except Exception as e:
-        return {"success": False, "error": str(e)}
+        return {"success": False, "result": None, "error": str(e)}
 
 
 def create_app() -> FastMCP: