"""Equation-based reservoir/neck sizing. SI flow, mm geometry, explicit input assumptions.
Not a coupled mushy-zone flow or calibrated feeding-distance solver.
"""
import math
import numpy as np
import trimesh
from filling import pipe_velocity

SOURCES=[
 'https://pmc.ncbi.nlm.nih.gov/articles/PMC9654928/',
 'https://archive.nptel.ac.in/content/storage2/courses/112107144/metalcasting/lecture13%2614.htm',
 'https://beckermann.lab.uiowa.edu/sites/beckermann.lab.uiowa.edu/files/2023-10/redbook.pdf',
]

def round_up(value):
    return math.ceil(value*10-1e-10)/10

def smallest_diameter(minimum,condition):
    lo=float(minimum);hi=lo
    if not 0<lo<=500:raise ValueError('Diameter must be within 0 to 500 mm.')
    while not condition(hi):
        if hi>=500:raise ValueError('No feasible feeder/channel diameter below 500 mm. Review input conditions.')
        hi=min(500.,hi*1.2)
    for _ in range(50):
        mid=(lo+hi)/2
        if condition(mid):hi=mid
        else:lo=mid
    result=round_up(hi)
    if not condition(result):result=round_up(result+.1)
    if result>500:raise ValueError('Required rounded diameter exceeds 500 mm.')
    return result

def cylinder_measures(diameter,height,neck_diameter=0,sections=64):
    # Match the actual polygonal STL, not the larger ideal circular volume.
    radius=diameter/2
    base=.5*sections*radius**2*math.sin(2*math.pi/sections)
    perimeter=2*sections*radius*math.sin(math.pi/sections)
    contact=.5*sections*(min(neck_diameter,diameter)/2)**2*math.sin(2*math.pi/sections)
    volume=base*height;cooling=perimeter*height+2*base-contact
    return volume,cooling,volume/cooling

def extend_to_pouring_plane(report,geometry,top):
    """A full-width open feeder reaches the same plane as the sprue in the exported mold.
    Extra height increases available volume and V/A; do not silently export a different feeder.
    """
    geometry['height_mm']=top-geometry['base_z_mm'];geometry['top_z_mm']=top
    vr,area,modulus=cylinder_measures(geometry['diameter_mm'],geometry['height_mm'],geometry['neck_diameter_mm'])
    report['minimum_height_mm']=report['height_mm'];report['height_mm']=geometry['height_mm']
    report['reservoir_volume_cm3']=vr/1000;report['reservoir_cooling_area_mm2']=area
    report['reservoir_modulus_mm']=modulus
    report['available_feed_after_reservoir_contraction_cm3']=(report['gross_usable_fraction_pct']-report['shrinkage_pct'])/100*vr/1000
    report['solidification_time_ratio']=report['riser_to_casting_constant_ratio']*(modulus/report['design_casting_modulus_mm'])**2
    report['checks']['volume']=report['available_feed_after_reservoir_contraction_cm3']>=report['required_feed_with_safety_cm3']
    report['checks']['solidification_time']=report['solidification_time_ratio']>=report['required_time_ratio']-1e-10
    report['all_sizing_constraints_met']=all(report['checks'].values())
    report['neck_capacity_is_conservative_at_minimum_reservoir_height']=True

def design_feeder(mesh,s,c,fs,gs,options,mask=None):
    beta=float(options.get('shrinkage_pct',6))/100
    eta=float(options.get('feeder_efficiency_pct',25))/100
    factor=float(options.get('feeding_safety',1.15))
    modulus_factor=float(options.get('modulus_factor',1.2))
    aspect=float(options.get('riser_height_ratio',1.5))
    feeding_time=float(options.get('feeding_time_s',120))
    peak=float(options.get('feeding_peak_factor',2))
    neck_min=float(options.get('neck_min_mm',3))
    neck_length_min=float(options.get('neck_length_mm',5))
    b_ratio=float(options.get('riser_constant_ratio',1))
    values=(beta,eta,factor,modulus_factor,aspect,feeding_time,peak,neck_min,neck_length_min,b_ratio)
    if not all(math.isfinite(v) for v in values):raise ValueError('All feeding inputs must be finite.')
    if not (0<beta<eta<=1 and factor>=1 and modulus_factor>1 and .5<=aspect<=3 and feeding_time>0 and peak>=1 and 0<neck_min<=500 and neck_length_min>0 and b_ratio>0):
        raise ValueError('Require 0 < shrinkage < feeder efficiency <= 100%, safety >= 1, modulus factor > 1, H/D 0.5–3, positive feeding time/neck/thermal constant ratio and peak factor >= 1.')
    volume=abs(float(mesh.volume));area=float(mesh.area)
    if min(volume,area)<=0:raise ValueError('A positive casting volume and cooling area are required.')
    global_modulus=volume/area
    local_modulus=0.
    if mask is not None and c.get('thick_point') is not None:
        # EDT selects the interior point, but the radius is measured to actual STL triangles.
        # Do not mistake cell-centre distance to an exterior voxel for exact surface distance.
        _,distance,_=trimesh.proximity.closest_point_naive(mesh,np.asarray([c['thick_point']],float))
        local_modulus=float(distance[0])/3
    casting_modulus=max(global_modulus,local_modulus)
    point=np.asarray(c['riser'],float)
    base_z=max(float(mesh.bounds[1,2])+neck_length_min,float(point[2])+neck_length_min)
    neck_length=base_z-float(point[2])
    # Lateral cooling only at the neck: both ends are metal contacts. Same B and n=2.
    polygon_modulus_per_d=math.cos(math.pi/64)/4
    neck=round_up(max(neck_min,casting_modulus/polygon_modulus_per_d))
    minimum=float(options.get('riser_min_mm',6))
    if not math.isfinite(minimum) or not 0<minimum<=500:raise ValueError('Invalid minimum reservoir diameter.')
    def neck_capacity(diameter,diameter_riser):
        height=aspect*diameter_riser
        # End-of-usable-withdrawal level; positive static head remains above neck entrance.
        head=(neck_length+height*(1-eta))/1000
        capillary=max(0.,-4*gs.surface_tension*math.cos(math.radians(gs.contact_angle_deg))/((diameter*math.cos(math.pi/64)/1000)))
        head-=capillary/(fs.liquid_density*9.80665)
        velocity=pipe_velocity(head,fs.liquid_density,fs.viscosity,diameter*math.cos(math.pi/64)/1000,neck_length/1000,fs.minor_loss,fs.roughness_mm/1000)
        cross=.5*64*(diameter/2000)**2*math.sin(2*math.pi/64)
        return cross*velocity*1e6
    # Neck size affects its own contraction demand and the reservoir's contact/cooling area.
    for _ in range(50):
        neck_volume=cylinder_measures(neck,neck_length+neck*.4+.05)[0]
        required=factor*beta*(volume+neck_volume)
        target_modulus=casting_modulus*modulus_factor/b_ratio**.5
        def feasible(d):
            vr,ar,mr=cylinder_measures(d,aspect*d,neck)
            return d>=neck and (eta-beta)*vr>=required and mr>=target_modulus
        diameter=smallest_diameter(max(minimum,neck),feasible)
        peak_flow=peak*beta*(volume+neck_volume)/feeding_time/1000
        larger=smallest_diameter(neck,lambda d:neck_capacity(d,diameter)>=peak_flow)
        if larger<=neck+1e-9:break
        neck=larger
    else:raise ValueError('Feeder/neck constraints did not converge.')
    height=aspect*diameter;vr,ar,mr=cylinder_measures(diameter,height,neck)
    # Both ends of neck contact liquid metal; exclude both end areas exactly.
    neck_modulus=(neck/4)*math.cos(math.pi/64)
    ratio=b_ratio*(mr/casting_modulus)**2
    checks={'volume':bool((eta-beta)*vr>=required),'solidification_time':bool(ratio>=modulus_factor**2-1e-10),'neck_solidification':bool(neck_modulus>=casting_modulus-1e-9),'neck_flow':bool(neck_capacity(neck,diameter)>=peak_flow)}
    report={'method':'volume-modulus-neck-balance-v1','shrinkage_pct':beta*100,'gross_usable_fraction_pct':eta*100,'volume_safety':factor,'casting_volume_cm3':volume/1000,'neck_volume_cm3':neck_volume/1000,'shrinkage_demand_cm3':beta*(volume+neck_volume)/1000,'required_feed_with_safety_cm3':required/1000,'available_feed_after_reservoir_contraction_cm3':(eta-beta)*vr/1000,'reservoir_volume_cm3':vr/1000,'reservoir_cooling_area_mm2':ar,'casting_global_modulus_mm':global_modulus,'casting_local_sphere_bound_mm':local_modulus,'design_casting_modulus_mm':casting_modulus,'reservoir_modulus_mm':mr,'required_modulus_mm':target_modulus,'chvorinov_exponent':2,'riser_to_casting_constant_ratio':b_ratio,'solidification_time_ratio':ratio,'required_time_ratio':modulus_factor**2,'diameter_mm':diameter,'height_mm':height,'neck_diameter_mm':neck,'neck_length_mm':neck_length,'neck_modulus_mm':neck_modulus,'feeding_time_s':feeding_time,'peak_factor':peak,'required_peak_feed_cm3_s':peak_flow,'neck_capacity_cm3_s':neck_capacity(neck,diameter),'checks':checks,'all_sizing_constraints_met':all(checks.values()),'spatial_feeding_verified':False,'sources':SOURCES,'assumptions':['Entire casting plus conservatively counted full neck is assigned to one feeder.','Efficiency is gross withdrawable reservoir fraction; its own solidification contraction is subtracted separately.','Chvorinov exponent n=2; B_riser/B_casting is an explicit input; B_neck/B_casting=1.','Open reservoir: full top area counted as cooling, base excludes neck contact. STL uses 64-sided cylinders.','Local sphere modulus is a geometric lower bound, not a resolved last-to-freeze thermal modulus.','Neck pressure check is a peak demand check, not a transient mushy-zone Darcy solver.','Spatial feeding distances, isolated hot spots and porosity are not certified by these sizing inequalities.']}
    geometry={'diameter_mm':diameter,'height_mm':height,'neck_diameter_mm':neck,'neck_length_mm':neck_length,'base_z_mm':base_z,'top_z_mm':base_z+height,'sections':64}
    c['feeding_geometry']=geometry
    return report