#include <windows.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
from stdint import *
import vipermath
import t, c


W: t.CDefine = 900
H: t.CDefine = 540
MAX_STEPS: t.CDefine = 400
ESCAPE: t.CDefine = 60.0
DISK_IN: t.CDefine = 2.6      # 吸积盘内缘 (单位: 史瓦西半径 Rs=1)
DISK_OUT: t.CDefine = 11.0    # 吸积盘外缘
CAM_DIST: t.CDefine = 22.0
CAM_ELEV: t.CDefine = 9.0     # 相机仰角(度): 小角度=侧视/透镜效果最强(星际穿越). 想更俯视可调到35
FOVF: t.CDefine = 0.55        # tan(半视场)

# ---------------- 向量 ----------------
class V3:
    x: t.CDouble
    y: t.CDouble
    z: t.CDouble

    def __init__(self, x: t.CDefine, y: t.CDefine, z: t.CDefine):
        self.x = x
        self.y = y
        self.z = z

    def __add__(self, other: 'V3') -> 'V3':
        return V3(self.x + other.x, self.y + other.y, self.z + other.z)

    def __sub__(self, other: 'V3') -> 'V3':
        return V3(self.x - other.x, self.y - other.y, self.z - other.z)

    def __mul__(self, other: t.CDouble) -> 'V3':
        return V3(self.x * other, self.y * other, self.z * other)

    def __xor__(self, other: 'V3') -> t.CDouble: # dot
        return self.x * other.x + self.y * other.y + self.z * other.z
    
    def cross(self, other: 'V3') -> 'V3':
        return V3(self.y * other.z - self.z * other.y, self.z * other.x - self.x * other.z, self.x * other.y - self.y * other.x)

    def __len__(self) -> t.CDouble:
        return vipermath.sqrt(c.Deref(self) ^ c.Deref(self))

    def nrm(self) -> 'V3':
        l = len(c.Deref(self)) # self.__len__()
        return self / l if l > (1e-12) else self

    def lerp(self, other: 'V3', t: t.CDouble) -> 'V3':
        return self + (other - self) * t

def clmp(x: t.CDouble, a: t.CDouble, b: t.CDouble) -> t.CDouble:
    return x if x >= a and x <= b else a if x < a else b

def smoothstep(a: t.CDouble, b: t.CDouble, x: t.CDouble) -> t.CDouble:
    t = clmp((x - a) / (b - a), 0, 1)
    return t * t * (3 - 2 * t)



# static inline V3 v(double x,double y,double z){V3 r={x,y,z};return r;}
# static inline V3 add(V3 a,V3 b){return v(a.x+b.x,a.y+b.y,a.z+b.z);}
# static inline V3 sub(V3 a,V3 b){return v(a.x-b.x,a.y-b.y,a.z-b.z);}
# static inline V3 scl(V3 a,double s){return v(a.x*s,a.y*s,a.z*s);}
# static inline double dot(V3 a,V3 b){return a.x*b.x+a.y*b.y+a.z*b.z;}
# static inline V3 cross(V3 a,V3 b){return v(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x);}
# static inline double len(V3 a){return sqrt(dot(a,a));}
# static inline V3 nrm(V3 a){double l=len(a); return l>1e-12?scl(a,1.0/l):a;}
# static inline V3 lerp(V3 a,V3 b,double t){return add(scl(a,1-t),scl(b,t));}
# static inline double clmp(double x,double a,double b){return x<a?a:(x>b?b:x);}
# static inline double smoothstep(double a,double b,double x){double t=clmp((x-a)/(b-a),0,1);return t*t*(3-2*t);}

# ---------------- 帧缓冲 / 相机 ----------------
fb: UINT32PTR
CAM: V3
FWD: V3
RIGHT: V3
UP: V3
# static V3 CAM, FWD, RIGHT, UP;
gT: t.CDouble = 0.0
# static volatile double gT = 0.0;

def setupCamera():
    e: t.CDouble = CAM_ELEV * vipermath.U_M_PI / 180.0
    CAM = V3(CAM_DIST * vipermath.cos(e), CAM_DIST * vipermath.sin(e), 0)
    FWD = (V3(0, 0, 0) - CAM).nrm()
    RIGHT = (FWD.cross(V3(0,1,0))).nrm()
    UP = RIGHT.cross(FWD)

# ---------------- 哈希 / 星空 ----------------
def hash3(x: t.CDouble, y: t.CDouble, z: t.CDouble) -> t.CDouble:
    n: t.CDouble = vipermath.sin(x*127.1 + y * 311.7 + z * 74.7) * 43758.5453
    return n - vipermath.floor(n)

def stars(d: V3) -> V3:
    _c: V3 = V3(0, 0, 0)
    for i in range(3):
        s: t.CDouble = 80.0 + i * 90.0
        px: t.CDouble = d.x * s
        py: t.CDouble = d.y * s
        pz: t.CDouble = d.z * s
        ix: t.CDouble = vipermath.floor(px)
        iy: t.CDouble = vipermath.floor(py)
        iz: t.CDouble = vipermath.floor(pz)
        h:  t.CDouble = hash3(ix, iy, iz)
        if h > 0.991:
            fx: t.CDouble = px - ix - 0.5
            fy: t.CDouble = py - iy - 0.5
            d2: t.CDouble = fx * fx + fy * fy + fz * fz
            br: t.CDouble = (h-0.991) / 0.009 * vipermath.exp(-d2 * 9.0) * 1.4
            _t: t.CDouble = hash3(iy, iz, ix)
            tint: V3 = V3(0.7,0.8,1.0) if (_t < 0.3) else (V3(1.0,0.8,0.6) if (_t > 0.8) else V3(1,1,1));
            _c = add(_c, scl(tint, br))
    return _c


# ---------------- 吸积盘着色 ----------------
static void diskShade(V3 cp,double rr,double t,V3*emit,double*alpha){
    double tt=(rr-DISK_IN)/(DISK_OUT-DISK_IN);
    double edge = smoothstep(0,0.07,tt)*smoothstep(0,0.18,1.0-tt);
    double ang=atan2(cp.z,cp.x);
    double w = 1.0/pow(rr,1.5);                 # 开普勒角速度
    double ph = ang - t*0.0011*w*9.0;           # 随时间旋转
    double n  = 0.5+0.5*sin(ph*5.0 - rr*1.6);
    double n2 = 0.5+0.5*sin(ph*13.0 + rr*0.8 + 1.7);
    double dens = (0.35+0.65*n)*(0.55+0.45*n2);
    double bright = (0.45+1.7*pow(1.0-tt,1.8))*edge;

    # 相对论多普勒增亮 (一侧更亮 + 微蓝移)
    V3 vdir = nrm(v(-cp.z,0,cp.x));
    double beta = 0.46/sqrt(rr);
    V3 toCam = nrm(sub(CAM,cp));
    double mu = dot(vdir,toCam);
    double dop = 1.0/(1.0-beta*mu);
    double beam = pow(dop,3.0);

    V3 c1=v(1.0,0.92,0.68), c2=v(1.0,0.42,0.13);
    V3 base=lerp(c1,c2,tt);
    base.z += (dop-1.0)*0.35; base.x -= (dop-1.0)*0.05; # 蓝移微调

    double inten = bright*dens*beam*1.25;
    *emit = scl(base,inten);
    *alpha = clmp(edge*dens*0.92,0,1);
}

# ---------------- 光子测地线积分 ----------------
static V3 trace(V3 dir){
    V3 pos=CAM, vel=dir;
    V3 hvec=cross(pos,vel);
    double h2=dot(hvec,hvec);                   # 角动量守恒
    double T=1.0; V3 col=v(0,0,0);
    double t=gT;

    for(int i=0;i<MAX_STEPS;i++){
        double r2=dot(pos,pos), r=sqrt(r2);
        if(r<1.0){ break; }                     # 落入视界 -> 黑
        if(r>ESCAPE){ col=add(col,scl(stars(nrm(vel)),T)); break; }

        double dt=clmp(r*0.08,0.03,2.2);
        # RK4 (a = -1.5 h2 pos / r^5)
        ACC(P) scl(P, -1.5*h2/pow(dot(P,P),2.5))
        V3 k1p=vel,            k1v=ACC(pos);
        V3 k2p=add(vel,scl(k1v,dt*0.5)), k2v=ACC(add(pos,scl(k1p,dt*0.5)));
        V3 k3p=add(vel,scl(k2v,dt*0.5)), k3v=ACC(add(pos,scl(k2p,dt*0.5)));
        V3 k4p=add(vel,scl(k3v,dt)),     k4v=ACC(add(pos,scl(k3p,dt)));
        V3 npos=add(pos,scl(add(add(k1p,scl(k2p,2)),add(scl(k3p,2),k4p)),dt/6.0));
        V3 nvel=add(vel,scl(add(add(k1v,scl(k2v,2)),add(scl(k3v,2),k4v)),dt/6.0));

        # 吸积盘平面 y=0 穿越检测
        if(pos.y*npos.y<0){
            double f=pos.y/(pos.y-npos.y);
            V3 cp=lerp(pos,npos,f);
            double rr=sqrt(cp.x*cp.x+cp.z*cp.z);
            if(rr>DISK_IN && rr<DISK_OUT){
                V3 emit; double a;
                diskShade(cp,rr,t,&emit,&a);
                col=add(col,scl(emit,T));
                T*=(1.0-a);
                if(T<0.01) break;
            }
        }
        pos=npos; vel=nvel;
    }
    return col;
}

# ---------------- 着色一个像素 ----------------
static uint32_t shade(int px,int py){
    double aspect=(double)W/H;
    double sx=(2.0*(px+0.5)/W-1.0)*aspect*FOVF;
    double sy=(1.0-2.0*(py+0.5)/H)*FOVF;
    V3 dir=nrm(add(FWD,add(scl(RIGHT,sx),scl(UP,sy))));
    V3 c=trace(dir);
    # 曝光 + 色调映射 + gamma
    double ex=1.15;
    c=scl(c,ex);
    c.x=c.x/(1.0+c.x); c.y=c.y/(1.0+c.y); c.z=c.z/(1.0+c.z);
    c.x=pow(c.x,1.0/2.2); c.y=pow(c.y,1.0/2.2); c.z=pow(c.z,1.0/2.2);
    int R=(int)(clmp(c.x,0,1)*255), G=(int)(clmp(c.y,0,1)*255), B=(int)(clmp(c.z,0,1)*255);
    return (R<<16)|(G<<8)|B;
}

# ---------------- 多线程 ----------------
typedef struct{int y0,y1;} Job;
static DWORD WINAPI worker(LPVOID p){
    Job*j=(Job*)p;
    for(int y=j->y0;y<j->y1;y++)
        for(int x=0;x<W;x++)
            fb[y*W+x]=shade(x,y);
    return 0;
}
static int NT=4;
static void renderFrame(void){
    HANDLE th[64]; Job jb[64];
    int per=H/NT;
    for(int i=0;i<NT;i++){
        jb[i].y0=i*per; jb[i].y1=(i==NT-1)?H:(i+1)*per;
        th[i]=CreateThread(0,0,worker,&jb[i],0,0);
    }
    WaitForMultipleObjects(NT,th,TRUE,INFINITE);
    for(int i=0;i<NT;i++) CloseHandle(th[i]);
}

# ---------------- 窗口 ----------------
static LRESULT CALLBACK WndProc(HWND h,UINT m,WPARAM w,LPARAM l){
    if(m==WM_DESTROY){PostQuitMessage(0);return 0;}
    if(m==WM_KEYDOWN && w==VK_ESCAPE){PostQuitMessage(0);return 0;}
    return DefWindowProc(h,m,w,l);
}

int WINAPI WinMain(HINSTANCE hI,HINSTANCE hP,LPSTR cmd,int show){
    SYSTEM_INFO si; GetSystemInfo(&si);
    NT=(int)si.dwNumberOfProcessors; if(NT<1)NT=1; if(NT>64)NT=64;

    fb=(uint32_t*)malloc((size_t)W*H*4);
    bmi.bmiHeader.biSize=sizeof(BITMAPINFOHEADER);
    bmi.bmiHeader.biWidth=W; bmi.bmiHeader.biHeight=-H; # 顶向下
    bmi.bmiHeader.biPlanes=1; bmi.bmiHeader.biBitCount=32; bmi.bmiHeader.biCompression=BI_RGB;
    setupCamera();

    WNDCLASS wc={0}; wc.lpfnWndProc=WndProc; wc.hInstance=hI;
    wc.lpszClassName="bh"; wc.hCursor=LoadCursor(0,IDC_ARROW);
    RegisterClass(&wc);
    RECT rc={0,0,W,H}; AdjustWindowRect(&rc,WS_OVERLAPPEDWINDOW,FALSE);
    HWND hwnd=CreateWindow("bh","Gargantua  (ESC quit)",WS_OVERLAPPEDWINDOW,
        CW_USEDEFAULT,CW_USEDEFAULT,rc.right-rc.left,rc.bottom-rc.top,0,0,hI,0);
    ShowWindow(hwnd,show);

    LARGE_INTEGER fq,t0,t1; QueryPerformanceFrequency(&fq);
    DWORD start=GetTickCount();
    MSG msg; int run=1; double fps=0; DWORD ftc=0; int fcnt=0;
    while(run){
        while(PeekMessage(&msg,0,0,0,PM_REMOVE)){
            if(msg.message==WM_QUIT){run=0;break;}
            TranslateMessage(&msg); DispatchMessage(&msg);
        }
        if(!run)break;
        gT=(double)(GetTickCount()-start);     # <<< 传入的时间 t (毫秒)
        QueryPerformanceCounter(&t0);
        renderFrame();
        QueryPerformanceCounter(&t1);
        HDC dc=GetDC(hwnd);
        SetDIBitsToDevice(dc,0,0,W,H,0,0,0,H,fb,&bmi,DIB_RGB_COLORS);
        ReleaseDC(hwnd,dc);

        fcnt++; DWORD now=GetTickCount();
        if(now-ftc>500){ fps=fcnt*1000.0/(now-ftc); fcnt=0; ftc=now;
            char buf[128]; sprintf(buf,"Gargantua  FPS:%.1f  threads:%d",fps,NT);
            SetWindowText(hwnd,buf);
        }
    }
    free(fb);
    return 0;
}