"""Шесть характеристик и обратный режим: полевые оценки и полный журнал поиска."""
import sys,json,math,os
ROOT=globals().get('TVV_ROOT','')
if ROOT in sys.path:sys.path.remove(ROOT)
sys.path.insert(0,ROOT)
from run_record import start,finish,save_json,sha
from series_session import FieldSession
from regimes import electrical,search
names=['run_characteristics_v3.py','series_session.py','build_series.py','geometry.py','physics.py','regimes.py','numerics.py','run_record.py','machine.json','winding.json','rotor-bh.json','stator-bh.json','characteristic-plan-v3.json','model-scope.json']
family=globals().get('TVV_CHARACTERISTIC','occ')
resume=globals().get('TVV_RESUME_DIR','')
if resume:
    folder=resume
    with open(folder+'/receipt.json') as f:record=json.loads(f.read())
    if record['family']!=family or 'completed_evaluations' not in record:raise ValueError('Несовместимый запуск характеристики')
    for name,digest in record['inputs_and_programs'].items():
        if sha(ROOT+'/'+name)!=digest or sha(folder+'/snapshot/'+name)!=digest:raise ValueError('Изменена версия: '+name)
else:
    folder,record=start(ROOT,globals().get('TVV_OUTPUT_ROOT',ROOT+'/runs'),'characteristic',names,'2D active-section educational characteristics; Lend excluded, DC resistance at 15 C')
    record.update(family=family,completed_evaluations={},field_points={})
def load(name):
    with open(folder+'/snapshot/'+name+'.json') as f:return json.loads(f.read())
machine=load('machine');plan=load('characteristic-plan-v3');rs=machine['resistance']['stator_phase_ohm'];omega=2*math.pi*50
session=FieldSession(globals(),machine,load('winding')['rows'],{'rotor':load('rotor-bh'),'stator':load('stator-bh')},folder,plan['mesh'])
cache={};rows=[];evaluations=[]
for key,item in record['completed_evaluations'].items():
    for name,digest in item['files'].items():
        if sha(folder+'/'+name)!=digest:raise ValueError('Повреждена полевая оценка')
    cache[key]=item['value'];evaluations.append(item['value'])
evaluations.sort(key=lambda row:row['evaluation'])
def field(current,n=None):
    n=plan['angles'] if n is None else n;key=','.join('%.12g'%v for v in current)+'/'+str(n)
    if key in cache:return cache[key]
    index=len(evaluations);values=[]
    for k in range(n):
        stem='eval-%04d-a%03d'%(index,k)
        if stem in record['field_points']:
            for name,digest in record['field_points'][stem].items():
                if sha(folder+'/'+name)!=digest:raise ValueError('Повреждена промежуточная точка')
            with open(folder+'/'+stem+'.json') as f:row=json.loads(f.read())
            if row['currents_dqf_A']!=current or abs(row['theta_deg']-360*k/n)>1e-9:raise ValueError('Другие входы промежуточной точки')
        else:
            row=session.solve(360*k/n,current,stem);save_json(folder+'/'+stem+'.json',row)
            record['field_points'][stem]={stem+ext:sha(folder+'/'+stem+ext) for ext in ('.json','.fem','.ans')}
            save_json(folder+'/receipt.json',record)
        values.append(row)
    mean={key:sum(r[key] for r in values)/n for key in ('energy_J','coenergy_J')}
    mean.update(currents_A=current[:],angles=n,lambda_dqf_Wb=[sum(r['lambda_dqf_Wb'][j] for r in values)/n for j in range(3)],
                mean_torque_Nm=sum(r['force_torque'][2] for r in values)/n,evaluation=index)
    mean.update(electrical(current,mean['lambda_dqf_Wb'],rs))
    mean['torque_identity_scaled_error']=abs(mean['mean_torque_Nm']-mean['torque_dq_Nm'])/max(abs(mean['mean_torque_Nm']),10000.)
    if mean['torque_identity_scaled_error']>.02:raise ValueError('Средний момент не согласован с dq-мощностью')
    evaluations.append(mean);cache[key]=mean;save_json(folder+'/evaluations.json',evaluations)
    files={}
    for k in range(n):files.update(record['field_points']['eval-%04d-a%03d'%(index,k)])
    record['completed_evaluations'][key]={'value':mean,'files':files}
    record['evaluations']=len(evaluations);save_json(folder+'/receipt.json',record);print('TVV_CURVE_EVALUATION',family,index,current)
    return mean
def solve(initial,bounds,residual,current_from_unknown):
    result=search(lambda x:residual(field(current_from_unknown(x))),initial,bounds,plan['scaled_residual'])
    save_json(folder+'/latest-search.json',result)
    if not result['converged']:raise ValueError('Недостижимый режим/поиск: '+result['reason'])
    value=field(current_from_unknown(result['x']));value=dict(value,search=result)
    return value,result['x']
def append(value):rows.append(value);save_json(folder+'/results.json',rows)
try:
    if family=='occ':
        for excitation in plan['occ_If_A']:append(field([0.,0.,excitation]))
    elif family=='nominal' or family=='regulation':
        voltage=plan['U_line_V'];tanphi=math.tan(math.acos(plan['power_factor']));initial=[1500.,-2.1]
        amplitudes=[math.sqrt(2)*plan['P_nominal_W']/(math.sqrt(3)*voltage*plan['power_factor'])] if family=='nominal' else [math.sqrt(2)*i for i in plan['current_rms_A'] if i>0]
        for amplitude in amplitudes:
            current=lambda x:[amplitude*math.cos(x[1]),amplitude*math.sin(x[1]),x[0]]
            residual=lambda v:[(v['U_line_rms_V']-voltage)/voltage,(v['Q_generator_var']-v['P_generator_W']*tanphi)/max((v['P_generator_W']**2+v['Q_generator_var']**2)**.5,1e6)]
            value,initial=solve(initial,[(0.,6000.),(-math.pi+.001,-.001)],residual,current);append(value)
    elif family=='scc':
        initial=[-1000.,-1.]
        for excitation in plan['scc_If_A']:
            current=lambda x:[x[0],x[1],excitation]
            value,initial=solve(initial,[(-60000.,0.),(-6000.,6000.)],lambda v:[v['ud_V']/plan['short_voltage_scale_V'],v['uq_V']/plan['short_voltage_scale_V']],current);append(value)
    elif family=='external':
        excitation=plan['fixed_If_A'];tanphi=math.tan(math.acos(plan['power_factor']));initial=[-2.1]
        for rms in plan['external_current_rms_A']:
            if not rms:append(field([0.,0.,excitation]));continue
            amplitude=math.sqrt(2)*rms;current=lambda x:[amplitude*math.cos(x[0]),amplitude*math.sin(x[0]),excitation]
            value,initial=solve(initial,[(-math.pi+.001,-.001)],lambda v:[(v['Q_generator_var']-v['P_generator_W']*tanphi)/max((v['P_generator_W']**2+v['Q_generator_var']**2)**.5,1e6)],current);append(value)
    elif family=='angle':
        initial=[-5000.,-100.];peak=plan['U_line_V']*math.sqrt(2/3)
        for angle in plan['load_angle_deg']:
            ud=peak*math.sin(math.radians(angle));uq=peak*math.cos(math.radians(angle))
            value,initial=solve(initial,[(-60000.,30000.),(-60000.,30000.)],lambda v:[(v['ud_V']-ud)/peak,(v['uq_V']-uq)/peak],lambda x:[x[0],x[1],plan['fixed_If_A']]);append(dict(value,target_delta_deg=angle))
    elif family=='u-curve':
        initial=[-3000.,-7000.]
        for excitation in plan['u_If_A']:
            value,initial=solve(initial,[(-60000.,30000.),(-60000.,0.)],lambda v:[(v['U_line_rms_V']-plan['U_line_V'])/plan['U_line_V'],(v['P_generator_W']-plan['u_power_W'])/plan['P_nominal_W']],lambda x:[x[0],x[1],excitation]);append(value)
        # Add both midpoints around the computed minimum, without changing U or P.
        ranked=sorted(rows,key=lambda r:r['currents_A'][2]);index=min(range(len(ranked)),key=lambda i:ranked[i]['I_phase_rms_A'])
        if index==0 or index==len(ranked)-1:raise ValueError('Минимум U-кривой не охвачен диапазоном возбуждения')
        for neighbor in (index-1,index+1):
            excitation=(ranked[index]['currents_A'][2]+ranked[neighbor]['currents_A'][2])/2
            value,initial=solve(ranked[index]['currents_A'][:2],[(-60000.,30000.),(-60000.,0.)],lambda v:[(v['U_line_rms_V']-plan['U_line_V'])/plan['U_line_V'],(v['P_generator_W']-plan['u_power_W'])/plan['P_nominal_W']],lambda x:[x[0],x[1],excitation]);append(value)
    else:raise ValueError('Неизвестная характеристика')
    # Independent angular refinement of the middle curve point, not just a software check.
    reference=rows[len(rows)//2];refined=field(reference['currents_A'],2*plan['angles'])
    differences=[abs(a-b)/max(abs(b),1.) for a,b in zip(reference['lambda_dqf_Wb'],refined['lambda_dqf_Wb'])]
    if max(differences)>.01:raise ValueError('Угловое среднее характеристики не сошлось')
    if family=='external' and reference['I_phase_rms_A'] and abs(refined['power_factor']-plan['power_factor'])>.001:
        raise ValueError('Уточнённая серия не сохраняет коэффициент мощности')
    if family=='scc' and max(abs(refined['ud_V']),abs(refined['uq_V']))/plan['short_voltage_scale_V']>plan['scaled_residual']:
        raise ValueError('Уточнённая угловая серия КЗ не удовлетворяет U=0 в заданном абсолютном допуске')
    save_json(folder+'/angular-check.json',{'base':reference,'refined':refined,'flux_scaled_changes':differences})
    finish(folder,record,'passed_characteristic_points_and_selected_angular_check')
except Exception as error:record['error']=str(error);save_json(folder+'/results.json',rows);finish(folder,record,'failed');raise
print('TVV_CHARACTERISTIC_RESULT',folder,family,len(rows))