"""Полевые величины и независимые интегральные проверки; встроенный Python."""
import math
from geometry import make_geometry
from physics import abc_to_dq0, maxwell_circle, torque_dq


def collect(api,machine,winding,theta_deg,currents,folder,stem,coverage=True,contour_samples=1440):
    names=[p+str(b)+' — ветвь статора' for p in 'ABC' for b in (1,2)]+['F — возбуждение ротора']
    circuits={n:list(api['mo_getcircuitproperties'](n)) for n in names}
    phase=[(circuits[p+'1 — ветвь статора'][2]+circuits[p+'2 — ветвь статора'][2])/2 for p in 'ABC']
    dq=abc_to_dq0(phase,math.radians(theta_deg))
    independent={n:0. for n in names};area_error=0.
    labels=make_geometry(machine,winding,theta_deg)['labels']
    for label in labels:
        if not label['circuit']:continue
        api['mo_selectblock'](label['x'],label['y'])
        area=api['mo_blockintegral'](5)
        integral_a_dv=api['mo_blockintegral'](1)
        api['mo_clearblock']()
        side=machine['rotor'] if label['group']==11 else machine['stator']
        expected=side['conductor_width_mm']*side['conductor_height_mm']*1e-6
        area_error=max(area_error,abs(area-expected)/expected)
        # Type 1 already includes depth: (N/S)*integral(A dV).
        independent[label['circuit']]+=label['turns']/expected*integral_a_dv
    api['mo_groupselectblock']()
    energy=api['mo_blockintegral'](2);coenergy=api['mo_blockintegral'](17)
    api['mo_clearblock']()
    ilambda=sum(circuits[n][0]*circuits[n][2] for n in names)
    api['mo_groupselectblock'](10);api['mo_groupselectblock'](11)
    torque=api['mo_blockintegral'](22)
    api['mo_clearblock']();api['mo_smooth']('off')
    contours=[];g=machine['geometry'];rr=g['rotor_body_diameter_mm']/2
    with open(folder+'/'+stem+'-gap.csv','w') as f:
        f.write('radius_mm,angle_rad,Bx_T,By_T\n')
        for fraction in ((.3,.5,.7) if contour_samples else ()):
            radius=rr+fraction*g['air_gap_mm'];samples=[]
            for k in range(contour_samples):
                angle=2*math.pi*(k+.5)/contour_samples
                value=api['mo_getpointvalues'](radius*math.cos(angle),radius*math.sin(angle))
                samples.append((angle,value[1],value[2]))
                f.write('%.12g,%.12g,%.12g,%.12g\n'%(radius,angle,value[1],value[2]))
            force=maxwell_circle(samples,radius/1000,g['stack_length_mm']/1000)
            contours.append({'radius_mm':radius,'samples':contour_samples,'Fx_N':force[0],'Fy_N':force[1],'torque_Nm':force[2]})
    result={'theta_deg':theta_deg,'currents_dqf_A':currents,'circuits':circuits,
        'lambda_abc_Wb':phase,'lambda_dq0_Wb':dq,'lambda_f_Wb':circuits[names[-1]][2],
        'independent_circuit_flux_Wb':independent,'max_copper_area_relative_error':area_error,
        'energy_J':energy,'coenergy_J':coenergy,'sum_i_lambda_J':ilambda,
        'torque_weighted_stress_Nm':torque,'torque_dq_estimate_Nm':torque_dq(dq,currents[:2]),
        'torque_dq_warning':'instantaneous estimate is not exact torque when spatial harmonics depend on rotor angle',
        'maxwell_contours':contours,'elements':api['mo_numelements'](),'nodes':api['mo_numnodes'](),
        'scope':'documented teaching model, not factory validation; no rated regime solved'}
    if coverage:
        _,centers,areas,groups=api['mo_getallelements']()
        induction=api['mo_getallb']();stats={}
        with open(folder+'/'+stem+'-field.csv','w') as f:
            f.write('x_mm,y_mm,group,Bx_T,By_T\n')
            for k in range(len(areas)):
                group=int(groups[k]);x,y=float(centers[k,0]),float(centers[k,1])
                bx,by=float(induction[k,0]),float(induction[k,1]);b=math.sqrt(bx*bx+by*by)
                if k%12==0 and x*x+y*y<(g['stator_back_diameter_mm']*.55)**2:
                    f.write('%.9g,%.9g,%d,%.9g,%.9g\n'%(x,y,group,bx,by))
                tag='stator' if group==20 else 'rotor' if group==10 and x*x+y*y>(g['rotor_inner_diameter_mm']/2)**2 else None
                if not tag:continue
                s=stats.setdefault(tag,{'area_mm2':0.,'low_area_mm2':0.,'high_area_mm2':0.,'max_bulk_B_T':0.,'energy_J':0.,'low_energy_J':0.,'high_energy_J':0.})
                values=api['mo_getpointvalues'](x,y)
                # Filled isotropic material: distinguish mean B from steel B.
                h=math.sqrt(values[5]**2+values[6]**2)
                fill=machine['stator']['stack_factor'] if tag=='stator' else 1
                steel=(b-(1-fill)*4e-7*math.pi*h)/fill
                low,high=machine['_bh_limits'][tag]
                area=float(areas[k]);e=values[4]*area*1e-6*g['stack_length_mm']/1000
                s['area_mm2']+=area;s['energy_J']+=e;s['max_bulk_B_T']=max(s['max_bulk_B_T'],b)
                if steel<low:s['low_area_mm2']+=area;s['low_energy_J']+=e
                if steel>high:s['high_area_mm2']+=area;s['high_energy_J']+=e
        for s in stats.values():
            for kind in ('low','high'):
                s[kind+'_area_fraction']=s[kind+'_area_mm2']/s['area_mm2']
                s[kind+'_energy_fraction']=s[kind+'_energy_J']/max(s['energy_J'],1e-12)
        result['material_coverage']=stats
        result['material_coverage_note']='steel B inferred under in-plane scalar mixture; not a measured law or general anisotropy model'
    return result