#calculate average field for some latitudes from numpy import genfromtxt,sqrt,mean,sin,pi from scipy.stats import mode import sys #data=genfromtxt('v1358ori_map_2013.txt',skip_header=4) #data=genfromtxt('LQHya_map_2010.dat') data=genfromtxt('BECet_map_CDS.dat') lim=60 #limiting latitude, calculate avg field above this lat=data[:,0] lon=data[:,1] Br=data[:,2] Bm=data[:,3] Ba=data[:,4] #Brig=data[:,5] ind=[] #indices of the grid points to be included will go here for i in range(0,len(lat)): if (lat[i] >= lim): ind.append(i) #weight the data points with the size of the grid point longitudes=[4,10,17,23,29,34,40,45,50,55,59,63,66,69,72,74,75,76,77,77,76,75,74,72,69,66,63,59,55,50,45,40,34,29,23,17,10,4] #number of longitudes on each latitude n=1876 nlat=38 weights=[] #w = sa * n_gridpoints / full_sa for i in range(0,len(Br)): if (abs(lat[i]) == 87.632): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) #sa=(sin(pi/2)-sin((pi/2-87.632*pi/180)*2))*(2*pi/4) sa=(sin(pi/2*(1-0/19))-sin(pi/2*(1-1/19)))*(2*pi/4) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 82.895): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-1/19))-sin(pi/2*(1-2/19)))*(2*pi/10) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 78.158): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-2/19))-sin(pi/2*(1-3/19)))*(2*pi/17) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 73.421): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-3/19))-sin(pi/2*(1-4/19)))*(2*pi/23) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 68.684): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-4/19))-sin(pi/2*(1-5/19)))*(2*pi/29) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 63.947): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-5/19))-sin(pi/2*(1-6/19)))*(2*pi/34) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 59.211): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-6/19))-sin(pi/2*(1-7/19)))*(2*pi/40) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 54.474): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-7/19))-sin(pi/2*(1-8/19)))*(2*pi/45) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 49.737): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-8/19))-sin(pi/2*(1-9/19)))*(2*pi/50) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 45.0): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-9/19))-sin(pi/2*(1-10/19)))*(2*pi/55) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 40.263): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-10/19))-sin(pi/2*(1-11/19)))*(2*pi/59) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 35.526): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-11/19))-sin(pi/2*(1-12/19)))*(2*pi/63) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 30.789): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-12/19))-sin(pi/2*(1-13/19)))*(2*pi/66) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 26.053): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-13/19))-sin(pi/2*(1-14/19)))*(2*pi/69) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 21.316): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-14/19))-sin(pi/2*(1-15/19)))*(2*pi/72) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 16.579): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-15/19))-sin(pi/2*(1-16/19)))*(2*pi/74) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 11.842): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-16/19))-sin(pi/2*(1-17/19)))*(2*pi/75) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 7.105): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-17/19))-sin(pi/2*(1-18/19)))*(2*pi/76) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) if (abs(lat[i]) == 2.368): #solid angle = (sin(latN)-sin(latS))*(lonE-lonW) sa=(sin(pi/2*(1-18/19))-sin(pi/2*(1-19/19)))*(2*pi/77) w=sa*n/(4*pi) #print(lat[i],w) weights.append(w) #sys.exit() #calculate rms (root mean square) of Br s=[] #sum of the squares Btot=[] #calculate also |B| Btotnw=[] Btotsq=[] #squared for rms wsum=[] #weighted sum of the polar field for i in range(0,len(Br)): s.append(weights[i]*Br[i]**2) Btot.append(weights[i]*sqrt(Br[i]**2+Bm[i]**2+Ba[i]**2)) #remove the sqrt?? Btotnw.append(sqrt(Br[i]**2+Bm[i]**2+Ba[i]**2)) #no weighting (for mode) Btotsq.append(weights[i]*(Br[i]**2+Bm[i]**2+Ba[i]**2)) #since its squared?? wsum.append(weights[i]*Br[i]) rms=sqrt(sum(s)/len(s)) rmstot=(sqrt(sum(Btotsq)/len(Btotsq))) modetot=float(mode(Btotnw)[0]) #weighted or not? moderad=float(mode(Br)[0]) #moderad=float(mode(wsum)[0]) Brad=[] #add here the field of the grid points to be included Bmer=[] Bazi=[] #Btot=[] for i in ind: #Brad.append(Br[i]) Brad.append(wsum[i]) Bmer.append(Bm[i]) Bazi.append(Ba[i]) #Btot.append(sqrt(Br[i]**2+Bm[i]**2+Ba[i]**2)) print('rms Brad: '+str(rms)+' kG') print('rms Btot: '+str(rmstot)+' kG') print('avg Brad: '+str(mean(Br))+' kG') print('avg |B|: '+str(mean(Btot))+' kG') print('mode Brad: '+str(moderad)+' kG') print('mode Btot: '+str(modetot)+' kG') print('avg Brad, lat > '+str(lim)+': '+str(mean(Brad))+' kG')