Calculating circle contact point and creating the 3-point arc in the vicinity of contact point

Question

Hii I am trying to calculate the contact point of two mating circles and for the calculated contact point I am trying to create 3 point arc in the vicinity of the contact point (for both the operation I am using contact point algorithm ). But the error in the second operation as i am storing those array in a text file and again importing that text file for the second operation. Is any better way to do that

Any help is appreciated

Thanks

    import numpy as np  
    from math import sqrt
    import math as math 
    fd =open('circle_input.txt','r')
    d=np.loadtxt(fd,delimiter=',',dtype={'names':('co1','col2','col3'),'formats':('float','float','float')})
    fd_1 =open('circle_input.txt','r')
    fd_2= open('circle_input.txt','r')
    d_1=np.loadtxt(fd_1,delimiter=',',dtype={'names':('co1','col2','col3'),'formats':('float','float','float')})
    d_2=np.loadtxt(fd_2,delimiter=',',dtype={'names':('co1','col2','col3'),'formats':('float','float','float')})
    temp1=d_1
    temp2=d_2

    for j in range (0,len(d)):

            if j==0:
                temp1[0]=temp1[j]
                temp1[j]=temp2[0]
                print(temp1)
            else:
                temp1[0]=temp1[j]
                temp1[j]=temp2[j-1]
                print(temp1)
                #print('\n')
            theta=0
            x0=d[j][0]+((d[j][2])*math.cos(theta))
            y0=d[j][1]+((d[j][2])*math.sin(theta))
            with open('output.txt','w') as f:
                for q in range(j+1,len(d)):
            if j!=q:
                x0,y0,r0,x1,y1,r1=d[j][0],d[j][1],d[j][2],d[q][0],d[q][1],d[q][2]
                dis=sqrt((x1-x0)**2+(y1-y0)**2)
                if dis<=(r0+r1):
                    a=(r0**2-r1**2+dis**2)/(2*dis)
                    h=sqrt(r0**2-a**2)
                    x2=x0+(a*(x1-x0)/dis)
                    y2=y0+(a*(y1-y0)/dis)
                    if dis==(r0+r1):
                        yc=y2-(h*(x1-x0)/dis)
                        xc=x2+(h*(y1-y0)/dis)
                        rc=((r0+r1)/4) # IMP Trick (#)
                        f.writelines(str(xc)+','+str(yc)+','+str(rc)+'\n') # important
                        f.writelines(str(d[j][0])+','+str(d[j][1])+','+str(d[j][2])+'\n')

                    else:# This part is passive
                        xc1=x2+(h*(y1-y0)/dis)
                        yc1=y2-(h*(x1-x0)/dis)
                        xc2=x2-(h*(y1-y0)/dis)
                        yc2=y2+(h*(x1-x0)/dis)
                                    print(str(xc)+','+str(yc)+','+str(rc)+'\n')
                                    print(str(d[j][0])+','+str(d[j][1])+','+str(d[j][2])+'\n')


            fd1 =open('output.txt','r')
            d1=np.loadtxt(fd1,delimiter=',',dtype={'names':('co1','col2','col3'),'formats':('float','float','float')})
              #command for upgrading thr oupput_contact.txt
            with open("output.txt") as f:
                with open("output_full.txt", "w") as f1:
                for line in f:
                    f1.write(line)
            open('output.txt', 'w')
            #
            with open('output_contact.txt','w') as f1:
                x0=np.zeros(len(d1))
                y0=np.zeros(len(d1)) 
                r0=np.zeros(len(d1)) 
                x1=np.zeros(len(d1))
                y1=np.zeros(len(d1)) 
                r1=np.zeros(len(d1)) 
                for p in range(0,len(d1)/2):
                    p=p+p
                    for q in range(p+1,p+2):
                    if p!=q:
                                     x0,y0,r0,x1,y1,r1=d1[p][0],d1[p][1],d1[p][2],d1[q][0],d1[q][1],d1[q][2]
                                     dis=sqrt((x1-x0)**2+(y1-y0)**2)                        
                                     if dis==(r0+r1):
                                         a=(r0**2-r1**2+dis**2)/(2*dis)
                                         h=sqrt(r0**2-a**2)
                                         x2=x0+(a*(x1-x0)/dis)
                                         y2=y0+(a*(y1-y0)/dis)
                                         if dis < (r0+r1):
                                             xc1=x2+(h*(y1-y0)/dis)
                                             yc1=y2-(h*(x1-x0)/dis)
                                             xc2=x2-(h*(y1-y0)/dis)
                                             yc2=y2+(h*(x1-x0)/dis)
                                             xc3=((x0+x1)/2)
                                             yc3=((y0+y1)/2)
                #f1.writelines(str(xc1)+','+str(yc1)+'\n'+str(xc2)+','+str(yc2)+'\n')
                print(str(xc1)+','+str(yc1)+','+str(xc2)+','+str(yc2)+','+str(xc3)+','+str(yc3)+'\n')




# circle_input.txt 
0,0,5
10,0,5
0,10,5