from numpy import *
from math import *
b1 = b2 = b3 = b4 = 1
ex = matrix('75,33,75;68,15,68;39,9,39;16,25,16;58,48,58;53,5,53;61,63,61;47,72,47;99,29,99;33,17,33;97,80,97;29,61,29;16,23,16;13,32,13;72,77,72;43,67,43;84,34,84;100,13,100;81,13,81;63,11,63')
Y = [0.00077385, 0.00042372, 8.8133, 4.5851, 9,4883, 0.0011336, 0.0011336, 0.10767, 0.43098, 0.80050, 0.00061111, 0.0031792, 44.0359, 0.01448, 0.0075917, 0.00026723, 0.000036466, 0.000095717, 0.000047435, 24.336]
m=20
n=4
eps=0.00001
X = zeros([m,n])
E = zeros([m,1])
B = zeros([n,1])
A = zeros([n,n])
bb1=bb2=bb3=bb4=k=0
while not ((abs((b1-bb1)/b1)<eps) and (abs((b2-bb2)/b2)<eps) and (abs((b3-bb3)/b3)<eps) and (abs((b4-bb4)/b4)<eps)):
    k=k+1
    print 'iteracya nomer', k
    print ''
    for element in range(0,m):
     n1 = sin(b4*ex[element,2])*((b1+b2*ex[element,0]**2+cos(b3*ex[element,1])))**(sin(b4*ex[element,2])-1)
     n2 = (sin(b4*ex[element,2])*((b1+b2*ex[element,0]**2+cos(b3*ex[element,1])))**(sin(b4*ex[element,2])-1))*ex[element,0]**2
     n3 = (sin(b4*ex[element,2])*((b1+b2*ex[element,0]**2+cos(b3*ex[element,1])))**(sin(b4*ex[element,2])-1))*(-sin(b3*ex[element,2])*ex[element,2])
     n4 = (sin(b4*ex[element,2])*((b1+b2*ex[element,0]**2+cos(b3*ex[element,1])))**(sin(b4*ex[element,2])))*log(b1+b2*ex[element,0]**2+cos(b3*ex[element,1]))*cos(ex[element,2]*b4)*ex[element,2]
     X[element,0] = n1
     X[element,1] = n2
     X[element,2] = n3
     X[element,3] = n4
     E[element,0] = Y[element]-(b1+b2*ex[element,0]**2+cos(b3*ex[element,1]))**sin(b4*ex[element,2])

    print 'matrix X:'
    print X
    A = dot(X.T,X)
    A = linalg.inv(A)
    Z = dot(X.T,E)
    B = dot(A,Z)
    print 'matrix B:'
    print B
    bb1=b1
    bb2=b2
    bb3=b3
    bb4=b4
    b1=b1+B[0]
    b2=b2+B[1]
    b3=b3+B[2]
    b4=b4+B[3]
    print b1,b2,b3,b4,k