#!/usr/bin/env python3
"""ABLAQ Studios — Radiative Environments, release 1.0.0.
Institutional research software. Scientific lead: Osama Khaled.
Python 3.10+, standard library only. No observations or weather files used.
Run: python3 -B radiative-study.v1.py --output results
License: CC0-1.0. Methods and references: https://ablaqstudios.com/data
"""
import argparse, csv, hashlib, json, math, pathlib, sys
from datetime import date
from collections import defaultdict
from time import perf_counter

VERSION = "1.0.0"
PI = math.pi
RATIOS = [1, 1.5, 2, 3, 4]
HEIGHTS = [2.5, 5, 7.5]
BEARINGS = list(range(0, 91, 15))
LATITUDES = [29, 32, 35]
MONTH_LENGTHS = [31,28,31,30,31,30,31,31,30,31,30,31]

def geometry(ratio, height):
    return math.sqrt(100/ratio), math.sqrt(100*ratio), height

def declination(day):
    # NOAA fractional-year approximation, evaluated at noon for each non-leap day.
    g = 2*PI*(day-1)/365
    return (.006918-.399912*math.cos(g)+.070257*math.sin(g)
            -.006758*math.cos(2*g)+.000907*math.sin(2*g)
            -.002697*math.cos(3*g)+.00148*math.sin(3*g))

def sun(latitude, decl, hour_angle):
    p = math.radians(latitude)
    e = -math.cos(decl)*math.sin(hour_angle)
    n = math.cos(p)*math.sin(decl)-math.sin(p)*math.cos(decl)*math.cos(hour_angle)
    z = math.sin(p)*math.sin(decl)+math.cos(p)*math.cos(decl)*math.cos(hour_angle)
    return e,n,z

def transmission(w,l,h,e,n,z,bearing):
    if z <= 0: return 0.
    b = math.radians(bearing)
    dx=h*(e*math.cos(b)-n*math.sin(b))/z
    dy=h*(n*math.cos(b)+e*math.sin(b))/z
    return max(0.,1-abs(dx)/w)*max(0.,1-abs(dy)/l)

def sky_exact(w,l,h):
    x,y=w/h,l/h
    return 2/(PI*x*y)*(
        .5*math.log((1+x*x)*(1+y*y)/(1+x*x+y*y))
        +x*math.sqrt(1+y*y)*math.atan(x/math.sqrt(1+y*y))
        +y*math.sqrt(1+x*x)*math.atan(y/math.sqrt(1+x*x))
        -x*math.atan(x)-y*math.atan(y))

def sky_angular(w,l,h,n=240):
    # Isotropic sky, cosine-weighted hemispherical integral. u=sin(altitude)^2.
    # Symmetry permits one azimuth quadrant. Midpoint product quadrature.
    total=0.
    angles=[(math.cos((i+.5)*PI/(2*n)),math.sin((i+.5)*PI/(2*n))) for i in range(n)]
    for j in range(n):
        u=(j+.5)/n
        distance=h*math.sqrt((1-u)/u)
        total+=sum(max(0,1-distance*c/w)*max(0,1-distance*s/l) for c,s in angles)
    return total/(n*n)

def sky_kernel(w,l,h,n=240):
    # Floor-to-opening area kernel, reduced using differences of coordinates.
    total=0.
    for i in range(n):
        x=(i+.5)*w/n
        for j in range(n):
            y=(j+.5)*l/n
            total+=(w-x)*(l-y)*h*h/(h*h+x*x+y*y)**2
    return total*4/(PI*n*n)

def solar_samples(latitude, step_minutes=10):
    groups=[[] for _ in range(12)]
    day=0
    for month,days in enumerate(MONTH_LENGTHS):
        for _ in range(days):
            day+=1; d=declination(day)
            for i in range(round(1440/step_minutes)):
                hour=(i+.5)*step_minutes/60
                e,n,z=sun(latitude,d,math.radians(15*(hour-12)))
                if z>0: groups[month].append((e,n,z,step_minutes/60))
    return groups

def integrate(groups,w,l,h,b):
    cb,sb=math.cos(math.radians(b)),math.sin(math.radians(b))
    rows=[]
    for samples in groups:
        beam=openbeam=hours=0.
        for e,n,z,dt in samples:
            dx=h*(e*cb-n*sb)/z;dy=h*(n*cb+e*sb)/z
            f=max(0.,1-abs(dx)/w)*max(0.,1-abs(dy)/l)
            beam+=f*z*dt;openbeam+=z*dt;hours+=f*dt
        rows.append((beam/openbeam,beam,openbeam,hours))
    return rows

def write_csv(path,rows):
    with path.open("w",newline="",encoding="utf-8") as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]),lineterminator="\n");w.writeheader()
        for row in rows: w.writerow({k:format(v,'.10g') if isinstance(v,float) else v for k,v in row.items()})

def svg_plot(path,title,description,series,xlabel,ylabel,xticks,yticks):
    # Publication figure, independent SVG artifact; all values supplied by calculations.
    from html import escape
    W,H=960,560;left,right,top,bottom=86,36,72,110
    pw,ph=W-left-right,H-top-bottom
    xmin,xmax=xticks[0][0],xticks[-1][0];ymin,ymax=yticks[0][0],yticks[-1][0]
    xp=lambda x:left+(x-xmin)/(xmax-xmin)*pw
    yp=lambda y:top+ph-(y-ymin)/(ymax-ymin)*ph
    a=[f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {W} {H}" role="img" aria-labelledby="t d"><title id="t">{escape(title)}</title><desc id="d">{escape(description)}</desc><rect width="960" height="560" fill="#f4f3ed"/><g font-family="Arial,sans-serif" fill="#142f29"><text x="28" y="34" font-size="20">{escape(title)}</text>']
    for y,label in yticks:
        a+=[f'<path d="M{left} {yp(y):.3f}H{W-right}" stroke="#d1d9cf"/><text x="{left-14}" y="{yp(y)+5:.3f}" font-size="15" text-anchor="end">{escape(str(label))}</text>']
    for x,label in xticks:
        a+=[f'<text x="{xp(x):.3f}" y="{top+ph+28}" font-size="15" text-anchor="middle">{escape(str(label))}</text>']
    a+=[f'<text x="{left+pw/2}" y="{H-52}" font-size="16" text-anchor="middle">{escape(xlabel)}</text><text x="22" y="{top+ph/2}" font-size="16" transform="rotate(-90 22 {top+ph/2})" text-anchor="middle">{escape(ylabel)}</text>']
    for i,(name,color,points) in enumerate(series):
        a+=[f'<polyline points="{" ".join(f"{xp(x):.3f},{yp(y):.3f}" for x,y in points)}" fill="none" stroke="{color}" stroke-width="3.5" stroke-dasharray="{["none","10 4","3 4"][i%3]}"/>']
        lx=left+i*260
        a+=[f'<path d="M{lx} {H-20}h28" stroke="{color}" stroke-width="4"/><text x="{lx+38}" y="{H-15}" font-size="14">{escape(name)}</text>']
    a+=['</g></svg>'];path.write_text(''.join(a)+'\n')

def main(out):
    out.mkdir(parents=True,exist_ok=False)
    monthly=[];annual=[];checks=[];samples={}
    for lat in LATITUDES:
        samples[lat]=solar_samples(lat)
        for r in RATIOS:
            for h in HEIGHTS:
                w,l,_=geometry(r,h);sv=sky_exact(w,l,h)
                for b in BEARINGS:
                    rows=integrate(samples[lat],w,l,h,b)
                    common=dict(latitude_deg_n=lat,aspect_ratio=r,wall_height_m=h,long_axis_bearing_deg=b)
                    for m,(t,beam,ob,hrs) in enumerate(rows,1):
                        monthly.append(dict(**common,month=m,beam_transmission=t,transmitted_cosine_hours=beam,open_cosine_hours=ob,sunlit_area_hours=hrs))
                    total=sum(x[1] for x in rows)/sum(x[2] for x in rows)
                    winter=sum(rows[i][1] for i in (11,0,1))/sum(rows[i][2] for i in (11,0,1))
                    summer=sum(rows[i][1] for i in (5,6,7))/sum(rows[i][2] for i in (5,6,7))
                    annual.append(dict(**common,annual_beam=total,winter_beam=winter,summer_beam=summer,sky_view_factor=sv))
    print(json.dumps({"phase":"solar","annual_cases":len(annual),"monthly_rows":len(monthly)}),flush=True)
    sky=[];maxangular=maxkernel=0.
    for r in [1,1.5,2,3,4,6,8]:
        for h in [1,2.5,5,7.5,10]:
            w,l,_=geometry(r,h);v=sky_exact(w,l,h)
            sky.append(dict(aspect_ratio=r,wall_height_m=h,width_m=w,length_m=l,area_m2=100,area_averaged_sky_view=v))
    for r,h in [(1,1),(1,5),(4,2.5),(4,7.5),(8,10)]:
        w,l,_=geometry(r,h);v=sky_exact(w,l,h)
        for n in [60,120,240]:
            angular=sky_angular(w,l,h,n);kernel=sky_kernel(w,l,h,n)
            checks.append(dict(kind="sky_quadrature",aspect_ratio=r,height=h,resolution=n,exact=v,angular=angular,kernel=kernel,angular_error=abs(angular-v),kernel_error=abs(kernel-v)))
            if n==240:maxangular=max(maxangular,abs(angular-v));maxkernel=max(maxkernel,abs(kernel-v))
    # Convergence for every latitude, aspect extremes, and both orientation extremes at h=5.
    maxsolar=0.
    for lat in LATITUDES:
        fine=solar_samples(lat,5)
        for r in [1,4]:
            for b in [0,90]:
                w,l,h=geometry(r,5)
                coarse=integrate(samples[lat],w,l,h,b);ref=integrate(fine,w,l,h,b)
                err=max(abs(a[0]-c[0]) for a,c in zip(coarse,ref))
                maxsolar=max(maxsolar,err)
                checks.append(dict(kind="temporal_convergence",latitude=lat,aspect_ratio=r,bearing=b,max_monthly_absolute_difference=err))
    # Positive-horizontal solar integral has a closed day-level answer at fixed declination.
    integral_error=0.
    for lat in LATITUDES:
        p=math.radians(lat)
        for day in [1,80,172,266,355]:
            d=declination(day);a=math.sin(p)*math.sin(d);b=math.cos(p)*math.cos(d)
            ws=math.acos(-a/b)
            exact=24/PI*(a*ws+b*math.sin(ws))
            num=sum(max(0,sun(lat,d,math.radians(15*((i+.5)/6-12)))[2])/6 for i in range(144))
            integral_error=max(integral_error,abs(num-exact))
    assert maxangular<0.0002, maxangular
    assert maxkernel<0.0005, maxkernel
    assert maxsolar<0.001, maxsolar
    assert integral_error<0.002, integral_error
    assert all(0<=r[k]<=1 for r in annual for k in ("annual_beam","winter_beam","summer_beam","sky_view_factor"))
    assert abs(sky_exact(10,10,10)-0.19982489569838746)<1e-12
    assert abs(sky_exact(5,20,5)-sky_exact(20,5,5))<1e-12
    assert transmission(10,10,0,1,1,1,0)==1
    assert transmission(10,10,5,0,0,1,0)==1
    # Square orientation invariance under quarter turn.
    for lat in LATITUDES:
        a=[x for x in annual if x['latitude_deg_n']==lat and x['aspect_ratio']==1 and x['wall_height_m']==5]
        assert abs(a[0]['annual_beam']-a[-1]['annual_beam'])<1e-12
    mix=[]
    for a in annual:
        if a['latitude_deg_n']!=32: continue
        for j in range(11):
            d=j/10
            mix.append(dict(**{k:a[k] for k in ('aspect_ratio','wall_height_m','long_axis_bearing_deg')},diffuse_fraction=d,
               annual_transmission=(1-d)*a['annual_beam']+d*a['sky_view_factor'],
               winter_transmission=(1-d)*a['winter_beam']+d*a['sky_view_factor'],
               summer_transmission=(1-d)*a['summer_beam']+d*a['sky_view_factor']))
    get=lambda ratio,bearing: next(x for x in annual if x['latitude_deg_n']==32 and x['aspect_ratio']==ratio and x['wall_height_m']==5 and x['long_axis_bearing_deg']==bearing)
    square=get(1,0);longns=get(4,0);longew=get(4,90)
    # Enumerate the non-dominated winter-admission / summer-rejection set, with fixed h=5.
    frontiers={}
    for j in [0,5,10]:
        rows=[a for a in mix if a["wall_height_m"]==5 and a["diffuse_fraction"]==j/10]
        frontiers[str(j/10)]=[a for a in rows if not any((b["winter_transmission"]>=a["winter_transmission"]-1e-12 and b["summer_transmission"]<=a["summer_transmission"]+1e-12 and (b["winter_transmission"]>a["winter_transmission"]+1e-12 or b["summer_transmission"]<a["summer_transmission"]-1e-12)) for b in rows)]
    diff=longns["winter_beam"]-square["winter_beam"]
    crossover=diff/(diff-(longns["sky_view_factor"]-square["sky_view_factor"]))
    summary=dict(version=VERSION,solar_cases=len(annual),monthly_rows=len(monthly),sky_cases=len(sky),mixture_cases=len(mix),latitude=32,height=5,square=square,long_ns=longns,long_ew=longew,winter_square_ns_crossover=crossover,frontiers=frontiers,
        accuracy=dict(max_monthly_10_vs_5_min=maxsolar,max_sky_angular_error=maxangular,max_sky_kernel_error=maxkernel,max_open_cosine_hour_error=integral_error))
    write_csv(out/"seasonal-monthly.v1.csv",monthly);write_csv(out/"seasonal-annual.v1.csv",annual)
    write_csv(out/"sky-view.v1.csv",sky);write_csv(out/"direct-diffuse.v1.csv",mix)
    (out/"radiative-summary.v1.json").write_text(json.dumps(summary,indent=2)+"\n")
    (out/"radiative-verification.v1.json").write_text(json.dumps(dict(status="passed",version=VERSION,verification_scope="Internal numerical comparison; no field validation or external review",tolerances=dict(solar_absolute_fraction=.001,angular_absolute_fraction=.0002,kernel_absolute_fraction=.0005,open_integral_cosine_hours=.002),**summary['accuracy'],comparisons=checks),indent=2)+"\n")
    sel=lambda r,b:[(a["month"],100*a["beam_transmission"]) for a in monthly if a["latitude_deg_n"]==32 and a["aspect_ratio"]==r and a["wall_height_m"]==5 and a["long_axis_bearing_deg"]==b]
    colors=['#173e34','#a16920','#52778b']
    svg_plot(out/"seasonal-exposure.v1.svg","Seasonal direct-beam transmission","Monthly cosine-weighted transmission, 32 degrees north, 100 square metre floor, five metre walls.",
        [(n,c,sel(r,b)) for n,c,r,b in zip(['Square, 0°','4:1, north–south','4:1, east–west'],colors,[1,4,4],[0,0,90])],"Month","Transmission (%)",[(1,"Jan"),(3,"Mar"),(5,"May"),(7,"Jul"),(9,"Sep"),(12,"Dec")],[(i,str(i)) for i in range(0,81,20)])
    svg_plot(out/"sky-access.v1.svg","Geometry and diffuse sky access","Area-averaged isotropic sky view factor versus wall height at three aspect ratios.",
        [(f"{r}:1 aspect ratio",c,[(a["wall_height_m"],100*a["area_averaged_sky_view"]) for a in sky if a["aspect_ratio"]==r]) for r,c in zip([1,4,8],colors)],"Wall height (m)","Sky view factor (%)",[(1,"1"),(2.5,"2.5"),(5,"5"),(7.5,"7.5"),(10,"10")],[(i,str(i)) for i in range(0,101,20)])
    svg_plot(out/"diffuse-tradeoff.v1.svg","Winter transmission under a mixed sky","Direct–diffuse mixture, 32 degrees north, December–February, five metre walls. Diffuse fraction is a scenario input, not a weather estimate.",
        [(n,c,[(a['diffuse_fraction'],100*a['winter_transmission']) for a in mix if a['aspect_ratio']==r and a['wall_height_m']==5 and a['long_axis_bearing_deg']==b]) for n,c,r,b in zip(['Square, 0°','4:1, north–south','4:1, east–west'],colors,[1,4,4],[0,0,90])],"Diffuse fraction of unobstructed horizontal exposure","Winter transmission (%)",[(x/4,str(x/4)) for x in range(5)],[(i,str(i)) for i in range(0,81,20)])
    files={p.name:dict(bytes=p.stat().st_size,sha256=hashlib.sha256(p.read_bytes()).hexdigest()) for p in sorted(out.iterdir())}
    files[pathlib.Path(__file__).name]=dict(bytes=pathlib.Path(__file__).stat().st_size,sha256=hashlib.sha256(pathlib.Path(__file__).read_bytes()).hexdigest())
    (out/"radiative-manifest.v1.json").write_text(json.dumps(dict(version=VERSION,publisher="ABLAQ Studios",author="ABLAQ Studios",scientific_lead=dict(name="Osama Khaled Ahmed Abdel-Qader",role="Scientific Research Lead",responsibility="Scientific direction, research decisions, and publication accountability"),method="Idealized radiative geometry; constant direct normal input; isotropic diffuse sky; no reflections or meteorological observations",year_modelled="Generic 365-day solar year, not historical weather",files=files),indent=2)+"\n")
    print(json.dumps(summary),flush=True)

if __name__=="__main__":
    ap=argparse.ArgumentParser(description=__doc__)
    ap.add_argument("--output",type=pathlib.Path,default=pathlib.Path("radiative-results"))
    args=ap.parse_args();main(args.output)