import sys
import math
max=[1.5, 1.5, 1.5]
min=[-1.5, -1.5, -1.5]

grain=[.02, .02, .02]

eye=[0., 0., -3.]

#light = eye
light=[0., 0., -3.]


# -0.2,0.4,-0.4,-0.4
#c=[-0.125,-0.256,0.847,0.0895]
#c=[-0.39554, 0.58898, 0.0, 0.0]

#c=[-0.450,-0.447,0.181,0.306]

#c=[0.645, -0.5, -0.113, -0.05]
c=[-0.213,-0.0410,-0.563,-0.560]

zarr=[[0] * int((max[0]-min[0])/grain[0])] * 3
zarri=0 #init stage

maxsearch=10;

maxcolor=511
maxiter=12
specular=.2
specularexp=1;

def printf(format, *args): sys.stdout.write (format % args)

print "P2 ", int((max[0]-min[0])/grain[0])-1, int((max[1]-min[1])/grain[1])-1, maxcolor
printf("#c:%s; light:%s; eye:%s; min:%s; max:%s; grain: %s\n", c, light, eye, min, max, grain)

def quat():
    global zarr, zarri, min, max, grain, eye
    viewy=max[1]

    zarri=0
    processrow(viewy) # top row, zarri=0
    zarri += 1
    
    viewy -= grain[1]
    processrow(viewy) #2nd row, zarri=1
    
    zarri += 1 #3rd row gonna be done, in zarri=2
    while(viewy > min[1]):
        viewy -= grain[1] #y moves down, pgm
        processrow(viewy)
        #just processed a row and stored in zarri... color previous one
        colorrow(viewy+grain[1]); #row ABOVE
        if (zarri==2):
            zarri=0
        else:
            zarri += 1
        #print "ZARR: " + `zarr`+"\n"
            

def colorrow(y):
    global zarri, zarr, grain, eye, maxcolor, min, max
    
    printf("#colorrow(%f)\n", y)
    if zarri is 0:
        zi=2
    else:
        zi = zarri-1
    if (zi == 0):
        zabove=2
        zbelow=1
    elif (zi == 1):
        zabove=0
        zbelow=2
    else:
        zabove=1
        zbelow=0
        
    x=min[0]
    i=0
    while (x < max[0]-grain[0]):
        if (zarr[zi][i]==max[2]):
            printf("%d " , maxcolor)
        else:
            if (x == min[0]):
                avec=[0-grain[0], 0, zarr[zi][i] - zarr[zi][i+1] ]
                bvec=[0, 2*grain[1], zarr[zabove][i] - zarr[zbelow][i]]
            else:
                avec=[2*grain[0], 0, zarr[zi][i+1] - zarr[zi][i-1]]
                bvec=[0, 2*grain[1], zarr[zabove][i] - zarr[zbelow][i]]
            #print "#zi: "+`zi`+"\n#avec: "+`avec`+" bvec: "+`bvec`

            #viewvec=[x-eye[0], y-eye[1], zarr[zi][i] - eye[2]]
            #try with light[]
            #lambertian reflection
            lightvec=[x-light[0], y-light[1], zarr[zi][i] - light[2]]
            nvec=cross(avec, bvec)
            lightvec=mkunitvec(lightvec)
            nvec=mkunitvec(nvec)

            phong=0
            if phong is 1:
                viewvec=[x-eye[0],  y-eye[1], zarr[zi][i] - eye[2]]
                viewvec=mkunitvec(viewvec)

                ldotn=dot(lightvec, nvec)
                rvec=nvec
                rvec[0]=rvec[0]*2*ldotn - lightvec[0]
                rvec[1]=rvec[1]*2*ldotn - lightvec[1]
                rvec[2]=rvec[2]*2*ldotn - lightvec[2]
                
                rvec=mkunitvec(rvec)
                col=maxcolor*ldotn + specular*(dot(rvec, viewvec)**specularexp)
            #lambert#col=ldotn*maxcolor
            #print "nvec, viewvec, col: ", nvec, viewvec, col

            #phong: R = 2(L . N)N - L
            #Kd * (N dot L) + Ks * (R dot V)^n
            else:
                col=dot(lightvec,nvec)*maxcolor
            col=abs(col)
            printf("%d ", int(col))
        x+=grain[0]
        i+=1
    printf("\n")

def dot(a, b):
    return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])

#a x b = (a2b3 - a3b2)i + (a3b1 - a1b3)j + (a1b2 - a2b1)k 
def cross(a, b):
    return [a[1]*b[2] - a[2]*b[1], a[2]*b[0] - a[0]*b[2], a[0]*b[1]-a[1]*b[0]]


# for each y in the viewing plane, iterate thru x and find quat-points
# for them all.
def processrow(viewy):
    global min, max, grain, zarr, zarri
    printf("#processrow(%f); zarri: %d\n", viewy, zarri)
    viewx=min[0]
    pix=0
    while (viewx < max[0]):
        zarr[zarri][pix] = processvector_orig(viewx, viewy)
        viewx += grain[0]
        pix+=1

# will end when distance betwen top/bottom is ~tiny
def processvector(x, y):
    global eye, min, grain, max, maxsearch
    
    point=[x, y, min[2] , 0] #point=x,y,min.z
    view=[0, 0, 1]
    
    i=0
    top=max[2]
    bottom=min[2]
    jp=-1
    while (1):
        #print "i, jp, top, bottom, p[2]:", i, jp, top, bottom, point[2]
        i+=1
        jp=jpoint(point)
        if jp != 1: #not in set, for forwards towards +z
            bottom = point[2]
            point[2] = point[2] + (top-point[2])/2.0
        else: #in set, go backwards
            top=point[2]
            point[2] = point[2] - (point[2]-bottom)/2.0
        #if (i>maxsearch and jp==1):
        if ((top-bottom)<.001):
            break
    if (top==max[2]): #probably blank
        return top
    return point[2]


# take in x, y on the viewing plane and return the z where the
# quat is at
def processvector_orig(x, y):
    global eye, min, grain, max,c
    point=[x, y, min[2], c[3]] #point=x,y,min.z
    while (point[2] < max[2]): #while point.z < max.z
        j=jpoint(point)
        if j is 1: return point[2]
        point[2]+=grain[2]
    return max[2]



# z[n+1] = z[n]^2 + c, c=cx + cy + cz + cw
def jpoint(pt):
    x=pt[0];
    y=pt[1]
    z=pt[2]
    w=pt[3]
    global c
    i=0
    #z[0]=x,y,z,w
    z=[x,y,z,0]
    znew=[0]*4    
    while ( i < maxiter):
        i+=1
        ztwo=2*z[0]
        znew[0]=( z[0]*z[0] - z[1]*z[1] - z[2]*z[2] - z[3]*z[3] + c[0] )
        znew[1]=( ztwo*z[1] + c[1] )
        znew[2]=( ztwo*z[2] + c[2] )
        znew[3]=( ztwo*z[3] + c[3] )
        #print "#znew: ",znew
        z=znew
        zmag=fake_mag(z)
        if (zmag > 8):
            return 0
    return 1

def fake_mag(vec):
    sum=0
    for i in vec:
        sum+=i**2
    return sum

def magnitude(vec):
    sum=0
    for i in vec:
        sum+=i**2
    return math.sqrt(sum)

def mkunitvec (vec):
    d=magnitude(vec)
    for i in range(len(vec)):
        vec[i]/=d
    return vec


if __name__=='__main__':
    quat()