11import numpy as np
22import time
33
4- # FIXME: Check this in its entirety
54def PCA (number : int , mass : np .ndarray , r : np .ndarray , Df : float , kf : float , tolerance : float ) -> tuple [bool , np .ndarray ]:
65 PCA_ok = True
76 n1 ,m1 ,rg1 ,x_cm ,y_cm ,z_cm , X ,Y ,Z = First_two_monomers (r , mass , number , Df , kf )
@@ -82,8 +81,6 @@ def PCA(number: int, mass: np.ndarray, r: np.ndarray, Df: float, kf: float, tole
8281 candidates *= 0
8382
8483
85- print (f"{ k = } " )
86- print ("WOOOOOOOOOOOW" )
8784 x_cm = (x_cm * m1 + X [k - 1 ]* m2 )/ (m1 + m2 )
8885 y_cm = (y_cm * m1 + Y [k - 1 ]* m2 )/ (m1 + m2 )
8986 z_cm = (z_cm * m1 + Z [k - 1 ]* m2 )/ (m1 + m2 )
@@ -123,12 +120,10 @@ def PCA_subcluster(N: int, N_subcluster: int, R: np.ndarray, DF: float, kf: floa
123120
124121 PCA_OK , data_new = PCA (number ,mass ,radius ,DF ,kf ,tolerance )
125122
126- # print(f"{PCA_OK = }")
127123 if i == 0 :
128124 acum = number
129125 for ii in range (number + 1 ):
130126 data [ii ,:] = data_new [ii ,:]
131- # ??
132127 i_orden [0 ,0 :2 ] = np .array ([1 , acum ])
133128 i_orden [0 ,2 ] = acum
134129 else :
@@ -152,8 +147,6 @@ def First_two_monomers(R: np.ndarray,M: np.ndarray,N: int,DF: float,kf:float) ->
152147 phi = np .arccos (2 * v - 1 )
153148 theta = 1
154149 phi = 1
155- # print(f"{theta = }, {phi = }")
156- # exit()
157150
158151 X [1 ] = X [0 ] + (R [0 ]+ R [1 ])* np .cos (theta )* np .sin (phi )
159152 Y [1 ] = Y [0 ] + (R [0 ]+ R [1 ])* np .sin (theta )* np .sin (phi )
@@ -163,17 +156,11 @@ def First_two_monomers(R: np.ndarray,M: np.ndarray,N: int,DF: float,kf:float) ->
163156 n1 = 2
164157
165158 rg1 = (np .exp (np .sum (np .log (R [:2 ]))/ 2 ))* np .power (n1 / kf ,1 / DF )
166- print (f"{ rg1 = } " )
167159
168160 x_cm = (X [0 ]* M [0 ]+ X [1 ]* M [1 ])/ (M [0 ] + M [1 ])
169161 y_cm = (Y [0 ]* M [0 ]+ Y [1 ]* M [1 ])/ (M [0 ] + M [1 ])
170162 z_cm = (Z [0 ]* M [0 ]+ Z [1 ]* M [1 ])/ (M [0 ] + M [1 ])
171163
172- print (f"{ X [0 ] = } , { Y [0 ] = } , { Z [0 ] = } " )
173- print (f"{ x_cm = } " )
174- print (f"{ y_cm = } " )
175- print (f"{ z_cm = } " )
176- print (f"{ M = } " )
177164 return n1 ,m1 ,rg1 ,x_cm ,y_cm ,z_cm , X ,Y ,Z
178165
179166def gamma_calc (rg1 : float ,rg2 : float ,rg3 : float ,n1 : int ,n2 : int ,n3 : int ) -> tuple [bool ,float ]:
@@ -284,8 +271,6 @@ def sticking_process(x: float,y: float,z: float,r: float,r_k: float, x_cm: float
284271
285272 alpha = np .arccos ((np .power (r1 ,2 ) + np .power (distance ,2 ) - np .power (r2 ,2 ))/ (2 * r1 * distance ))
286273 r0 = r1 * np .sin (alpha )
287- # alpha_0 = acos((r1**2.+distanc**2.-r2**2.)/(2.*r1*distanc))
288- # r0 = r1*sin(alpha_0)
289274
290275 # AmBdC = (A+B)/C
291276 AmBdC = (a + b )/ c
@@ -306,7 +291,8 @@ def sticking_process(x: float,y: float,z: float,r: float,r_k: float, x_cm: float
306291 return x_k , y_k , z_k , r0 , x0 ,y0 ,z0 ,i_vec , j_vec
307292
308293def sticking_process2 (x0 , y0 , z0 , r0 ,i_vec ,j_vec ):
309- theta = 2 * np .pi * np .random .rand ()
294+ u = np .random .rand ()
295+ theta = 2 * np .pi * u
310296
311297 x_k = x0 + r0 * np .cos (theta )* i_vec [0 ]+ r0 * np .sin (theta )* j_vec [0 ]
312298 y_k = y0 + r0 * np .cos (theta )* i_vec [1 ]+ r0 * np .sin (theta )* j_vec [1 ]
0 commit comments