Changeset 4099
 Timestamp:
 Aug 18, 2019 6:06:56 PM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/G2shapes.py
r4098 r4099 800 800 return dr; 801 801 802 # Return nonzero distances one coordinate & the rest 803 804 def get_drs(xyz,XYZ): 805 806 return np.sqrt(np.sum((XYZxyz)**2,axis=1)) 807 802 808 # Find center of beads within a radii 803 809 … … 806 812 num_beads = len(aList_beads_x) 807 813 808 xsum = 0.0 809 ysum = 0.0 810 zsum = 0.0 811 count_beads = 0.0 812 813 i = 0 814 while i < num_beads: 814 # xsum = 0.0 815 # ysum = 0.0 816 # zsum = 0.0 817 # count_beads = 0.0 818 # 819 # i = 0 820 # while i < num_beads: 821 # 822 # dr = get_dr(x,y,z,aList_beads_x[i],aList_beads_y[i],aList_beads_z[i]) 823 # 824 # if dr > radii_1 and dr < radii_2: 825 # count_beads = count_beads + 1.0 826 # xsum = xsum + float(aList_beads_x[i]) 827 # ysum = ysum + float(aList_beads_y[i]) 828 # zsum = zsum + float(aList_beads_z[i]) 829 # 830 # i = i + 1 831 # 832 # if count_beads > 0.0: 833 # xsum = xsum/count_beads 834 # ysum = ysum/count_beads 835 # zsum = zsum/count_beads 836 # delta = (xsum  x)**2 + (ysum  y)**2 + (zsum  z)**2 837 # delta = math.sqrt(delta) 838 # else: 839 # delta = 0.0 815 840 816 dr = get_dr(x,y,z,aList_beads_x[i],aList_beads_y[i],aList_beads_z[i]) 817 818 if dr > radii_1 and dr < radii_2: 819 count_beads = count_beads + 1.0 820 xsum = xsum + float(aList_beads_x[i]) 821 ysum = ysum + float(aList_beads_y[i]) 822 zsum = zsum + float(aList_beads_z[i]) 823 824 i = i + 1 825 826 if count_beads > 0.0: 827 xsum = xsum/count_beads 828 ysum = ysum/count_beads 829 zsum = zsum/count_beads 830 delta = (xsum  x)**2 + (ysum  y)**2 + (zsum  z)**2 831 delta = math.sqrt(delta) 832 else: 833 delta = 0.0 834 841 XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T 842 xyz = np.array([x,y,z]) 843 drs = get_drs(xyz,XYZ) 844 sumXYZ = np.array([XYZ[i] for i in range(num_beads) if radii_1<drs[i]<radii_2]) 845 count_beads = sumXYZ.shape[0] 846 847 delta = 0.0 848 if count_beads: 849 delta = np.sqrt(np.sum(((np.sum(sumXYZ,axis=0)/count_beads)xyz)**2)) 850 835 851 return delta; 836 852 … … 1027 1043 num_hist = len(aList_r) 1028 1044 max_dr = (float(num_hist)1.0)*hist_grid 1029 num_beads = len(aList_beads_x) 1030 1031 i = 0 1032 while i < num_hist: 1033 aList_pr_model_test2[i] = 0.0 1034 i = i + 1 1035 1036 i = 0 1037 while i < num_beads: 1038 1039 if i != ii: 1040 x2 = float(aList_beads_x[i]) 1041 y2 = float(aList_beads_y[i]) 1042 z2 = float(aList_beads_z[i]) 1043 dr = get_dr(x1,y1,z1,x2,y2,z2) 1044 dr = min(dr,max_dr) 1045 dr_grid = dr/hist_grid 1046 int_dr_grid = int(dr_grid) 1047 int_dr_grid = max(int_dr_grid,0) 1048 int_dr_grid = min(int_dr_grid,num_hist2) 1049 ip_low = int_dr_grid 1050 ip_high = ip_low + 1 1051 ip_high_frac = dr_grid  float(int_dr_grid) 1052 ip_low_frac = 1.0  ip_high_frac 1053 aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac 1054 aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac 1055 1056 i = i + 1 1057 1045 # num_beads = len(aList_beads_x) 1046 1047 aList_pr_model_test2 = num_hist*[0.0,] 1048 1049 XYZ = np.array([aList_beads_x,aList_beads_y,aList_beads_z]).T 1050 xyz = np.array([x1,y1,z1]) 1051 drs = get_drs(xyz,XYZ) 1052 drs_grid = np.where(drs>max_dr,max_dr,drs)/hist_grid 1053 int_drs_grid = np.array(drs_grid,dtype=np.int) 1054 int_drs_grid = np.where(int_drs_grid>num_hist2,num_hist2,int_drs_grid) 1055 ip_lows = int_drs_grid 1056 ip_highs = ip_lows + 1 1057 ip_high_fracs = drs_grid  int_drs_grid 1058 ip_low_fracs = 1.0  ip_high_fracs 1059 for ip_low in ip_lows: 1060 aList_pr_model_test2[ip_low] += ip_low_fracs[ip_low] 1061 for ip_high in ip_highs: 1062 aList_pr_model_test2[ip_high] += ip_high_fracs[ip_high] 1063 1064 # i = 0 1065 # while i < num_hist: 1066 # aList_pr_model_test2[i] = 0.0 1067 # i = i + 1 1068 # 1069 # i = 0 1070 # while i < num_beads: 1071 # 1072 # if i != ii: 1073 # x2 = float(aList_beads_x[i]) 1074 # y2 = float(aList_beads_y[i]) 1075 # z2 = float(aList_beads_z[i]) 1076 # dr = get_dr(x1,y1,z1,x2,y2,z2) 1077 # dr = min(dr,max_dr) 1078 # dr_grid = dr/hist_grid 1079 # int_dr_grid = int(dr_grid) 1080 # int_dr_grid = max(int_dr_grid,0) 1081 # int_dr_grid = min(int_dr_grid,num_hist2) 1082 # ip_low = int_dr_grid 1083 # ip_high = ip_low + 1 1084 # ip_high_frac = dr_grid  float(int_dr_grid) 1085 # ip_low_frac = 1.0  ip_high_frac 1086 # aList_pr_model_test2[ip_low] = float(aList_pr_model_test2[ip_low]) + ip_low_frac 1087 # aList_pr_model_test2[ip_high] = float(aList_pr_model_test2[ip_high]) + ip_high_frac 1088 # 1089 # i = i + 1 1090 # 1058 1091 return; 1059 1092
Note: See TracChangeset
for help on using the changeset viewer.