Twins

# -*- coding: utf-8 -*-
"""
Created on Thurs Sept 07 2015

@author: Tom Babbitt
"""
import random,math,cmath
from PIL import Image

class Fd(object):
    def __init__(self,*args):
        if len(args) < 3:
            if type(args[0]) == complex:
                real = args[0].real
                imag = args[0].imag
                if len(args) == 2:
                    funk = args[1]
            else:
                real = args[0]
                imag = 0.
                funk = 0.
        else:
            real = args[0]
            imag = args[1]
            funk = args[2]

        self.real = float(real)
        self.imag = float(imag)
        self.funk = float(funk)

    def distance_real(self):
        if (-self.funk**2 - 2*self.imag*self.funk - self.imag**2) < 0.:
            return Fd(math.sqrt(abs(-self.funk**2 - 2*self.imag*self.funk - self.imag**2)),.33333,.33333)
        else:
            return Fd(math.sqrt(-self.funk**2 - 2*self.imag*self.funk - self.imag**2),0.,0.)

    def distance_imag(self):
        if (-self.funk**2 - 2*self.real*self.funk + self.real**2) < 0.:
            return Fd(0.,math.sqrt(abs(-self.funk**2 - 2*self.real*self.funk + self.real**2)),.33333)
        else:
            return Fd(0.,math.sqrt(-self.funk**2 - 2*self.real*self.funk + self.real**2),0.)

    def distance_funk(self):
        if (-self.imag**2 - 2*self.real*self.imag + self.real**2) < 0.:
            return Fd(0.,.33333,math.sqrt(abs(-self.imag**2 - 2*self.real*self.imag + self.real**2)))
        else:
            return Fd(0.,0.,math.sqrt(-self.imag**2 - 2*self.real*self.imag + self.real**2))

    def distance(self):
        a = self.distance_real()
        b = a.Add(a,self.distance_imag())
        c = b.Add(b,self.distance_funk())
        return c

    def Add(self,Fd1,Fd2):
        real = Fd1.real + Fd2.real
        imag = Fd1.imag + Fd2.imag
        funk = Fd1.funk + Fd2.funk
        return Fd(real,imag,funk)

    def Subtract(self,Fd1,Fd2):
        real = Fd1.real - Fd2.real
        imag = Fd1.imag - Fd2.imag
        funk = Fd1.funk - Fd2.funk
        return Fd(real,imag,funk)

    def Multiply(self,Fd1,Fd2):
        a,b,c = Fd1.real,Fd1.imag,Fd1.funk
        d,e,f = Fd2.real,Fd2.imag,Fd2.funk
        if (a,b,c) == (0.,0.,0.):
            return Fd(0.,e+d,f+d)
        elif (d,e,f) == (0.,0.,0.):
            return Fd(0.,a+b,a+c)

        real = (a*d) - (b*e) - (c*f)
        imag = (a*e) + (b*d)
        funk = (a*f) + (c*d)
        imunk= (b*f) + (c*e)
        if imunk < 0.:
            real = real - imunk
        else:
            real = real + imunk
        return Fd(real,imag,funk)

    def Conjigate(self,Fd1):
        real = Fd1.real
        imag = -Fd1.imag
        funk = -Fd1.funk
        return Fd(real,imag,funk)

    def Reciprocal(self,Fd1):
        Fdr = self.conjigate(Fd1)
        divisor = Fd1.real**2 - Fd1.imag**2 - Fd1.funk**2
        return Fd(Fdr.real/divisor , Fdr.imag/divisor , Fdr.funk/divisor)

    def Divide(self,Fd1,Fd2):
        if (Fd1.real, Fd1.imag, Fd1.funk) == (0.,0.,0.):
            return Fd(0.,Fd2.imag+(1./Fd2.real),Fd2.funk+(1./Fd2.real))
        elif (Fd2.real, Fd2.imag, Fd2.funk) == (0.,0.,0.):
            return Fd(0.,Fd1.imag+(1./Fd1.real),Fd1.funk+(1./Fd1.real))
        divis = Reciprocal(Fd2)
        return Multiply(Fd1,divis)

    def Power(self,Fd1,n_pow):
        Fd2 = self.Multiply(Fd1,Fd1)
        if n_pow == 2:
            return Fd2
        else:
            for i in range(n_pow-1):
                Fd2 = self.Multiply(Fd2,Fd1)
            return Fd2

    def Absolute(self,Fd1):
        return math.sgrt(Fd1.real**2 - Fd1.imag**2 - Fd1.funk**2)
def mandle(a1):
    j,c = a1[0],a1[1]
    j = j.Power(j,4)
    j = j.Add(j,c)
    return ([j,c])

class IFS(object):
    def __init__(self,Fn,Io):
        self.Fn = Fn
        self.Io = Io
        self.Inext = Io
    def __iter__(self):
        return self
    def next(self):
        if self.Fn:
            self.Inext = self.Fn(self.Inext)
        return self.Inext

def checkval(Valu):
    vrl,vimg,vfnk = Valu.real,Valu.imag,Valu.funk
    if math.isnan(vrl):
        Valu = Fd(1.,vimg,vfnk)
        vrl = 1.
    if math.isnan(vimg):
        Valu = Fd(vrl,1.,vfnk)
        vimg = 1.
    if math.isnan(vfnk):
        Valu = Fd(vrl,vimg,1.)
        vimg = 1.
    if vrl > 1.:
        Valu = Fd(1.,vimg,vfnk)
        vrl = 1.
    elif vrl < -1.:
        Valu = Fd(-1.,vimg,vfnk)
        vrl = -1.
    if vimg > 1.:
        Valu = Fd(vrl,1.,vfnk)
        vimg = 1.
    elif vimg < -1.:
        Valu = Fd(vrl,-1.,vfnk)
        vimg = -1.
    if vfnk > 1.:
        Valu = Fd(vrl,vimg,1.)
        vfnk = 1.
    elif vfnk < -1.:
        Valu = Fd(vrl,vimg,-1.)
        vfnk = -1.
    return Valu

superimg = [[0. for i in range(4098)] for j in range(4098)]
superimg2 = [[0. for i in range(4098)] for j in range(4098)]

def imgux(gray,fname):
    im = Image.new('L',(len(gray[int(len(gray)/2)]),len(gray)))
    #for c,i in enumerate(gray):
    #    for k,l in enumerate(i):
    #        #print(c,k,l,)
    #        im.putpixel((c,k),int(l))
    [[im.putpixel((c,k),int(l)) for k,l in enumerate(i)] for c,i in enumerate(gray)]
    im.save(fname)

for FD in (Fd(-.7,0.1333333333333333,.2),Fd(-.7,0.1666666666666667,.2),Fd(-.7,-0.1666666666666667,.2),Fd(-.7,-0.1333333333333333,.2)):
    #((34-30)/30.,35,25,26)#(10,11,12,13,14,15,16,17,18,19,20,21,26)   #0.6333333333333333,.4 ,.2
    for extraitir in range(60):
        rnum = Fd(random.uniform(-1.,1.),random.uniform(-1.,1.),random.uniform(-1.,1.))
        #rnum = Fd(((extraitir-30.)/ 30.),((extraitir-30.)/ 30.),((extraitir-30.)/ 30.))
        newIFS = IFS(mandle,[rnum,FD])
        for itercount in range(2048):
            nextnum = newIFS.next()
            nextnum[0] = checkval(nextnum[0])
            superimg[int(nextnum[0].imag * 2048.)+ 2047][int(nextnum[0].real * 2048.)+ 2047]  = 225#int((i+1) * 4.25)
            superimg2[int(nextnum[0].funk * 2048.)+ 2047][int(nextnum[0].real * 2048.)+ 2047] = 225#int((i+1) * 4.25)

imgux(superimg,'c:\\tom2\\twins\\realimag.png')
imgux(superimg2,'c:\\tom2\\twins\\realfunk.png')