---------------------------------------------------------------------- Ulysses Mission EPAC - Energetic Particle Composition Instrument Algorithm Descriptions ---------------------------------------------------------------------- INTRODUCTION: This bundle contains EPAC pulse height analysis data in 1 day (24 hr) files. Data files are provied in both binary and ascii formats which contain the same information but in a slightly different format. The file naming convention is EV92xxx (.BIN binary) or (.TAB ascii table) where the xxx is the day of 1992. This file contains BIN_DATA, CHANTOEN, and EPAC_EV IDL and translated Python code algorithms to provide a means to use the IDL code provided by the EPAC Instrument Team. BIN_DATA describes the contents of the binary data files. CHANTOEN will convert EPAC channel numbers to energy (MeV). EPAC_EV reads all available PHA data and generates a postscript plot file: "epac-ev.ps". -------------------------------------------------------------------------- BIN_DATA -------------------------------------------------------------------------- This file Originally Delievered in IDL EV data files ============= The file names for ULYSSES-KEP Experiment EV data files are: EVyearday.BIN example: EV92013.BIN Data are written in a 3 bytes per event data stream with a variable record length. dE=0 and dx=0 records are excluded. Each file contains the events of one day. The date of the corresponding day is part of the file name. Contents of the 3 bytes: ---------------------------------------------------------------- Byte One: E (channel number) Byte Two: dE (channel number) Byte Three: ER * 128 + SpID * 16+ TID * 4 + SeID with: ER = Energy Range. 0 or 1, for lower or upper part E=E+256*SpID [ idl code: E=E+255*SpID ] SpID = Species ID, [0,1,2,3,4] for out and HS 1,2,3,4 TID = Telescope ID, [0,1,2,3] for Telescope 1,2,3,4 SeID = Sector ID, [0,1,2,3] for Sector 1,2,3,4 ------------------------------------------------------------------ The data are written with the following IDL code: ------------------------------------------------------------------ barr=BYTARR(3,vcounts) barr(0)=byte(E(index)) barr(1)=byte(dE(index)) barr(2)=byte(ER(index)*128+SpID(index)*16+TID(index)*4+SeID(index)) with: index: array index of valid data points vcounts: number of valid data points File opening: under VMS: openw,1,outputfile,/NONE, /BLOCK under other systems: openw,1,outputfile Writing: WRITEU,1,barr -------------------------------------------------------------------- Notes: E and dE channel numbers can be converted to energy using the IDL code contained in CHANTOEN.TXT. Python translation import os # Declare Path input_BIN_directory = "" # Assign Constants ENERGY_RANGE_UPPER = 1 ENERGY_RANGE_LOWER = 0 # Convert data into 3 bytes format def convert_to_3_bytes(channel_number, species_id, telescope_id, sector_id, energy_range=ENERGY_RANGE_LOWER): byte_one = channel_number & 0xFF byte_two = (channel_number >> 8) & 0xFF byte_three = (energy_range << 7) + (species_id << 4) + (telescope_id << 2) + sector_id return bytes([byte_one, byte_two, byte_three]) # Read data from binary files def read_bin_files(directory): for filename in os.listdir(directory): if filename.endswith(".BIN"): with open(os.path.join(directory, filename), "rb") as file: data = file.read() # Output 3 bytes format to txt file output_filename = os.path.splitext(filename)[0] + ".txt" output_path = os.path.join(directory, output_filename) with open(output_path, "w") as output_file: output_file.write(data.decode('utf-8', errors='ignore')) read_bin_files(input_BIN_directory) -------------------------------------------------------------------------- CHANTOEN -------------------------------------------------------------------------- ; NAME: CHANTOEN ; ; PURPOSE: Channel to Energy conversion for EPAC EV-Matrix ; ; CATEGORY: Element Distr. & Energy calc. ; ; CALLING SEQUENCE: en=CHANTOEN(chan,xybit,hmodbit) ; ; INPUTS: chan y:0-255 (dE), x:0-511 (E) EV channel ; or array of these numbers ; xybit 0: Xaxis (E, B-Det) 1: Yaxis (dE, A-Det) ; hmodbit 0:IMOD 1:HMOD (value of column ER for E, 0 for dE) ; ; OUTPUTS: En FLT or FLTARR respective energies (in MeV) ; ; RESTRICTIONS: ; ; PROCEDURE: ; Factors from experiment description, p.4-26a ; (recalibrated values, oct. 91) ; ; MODIFICATION HISTORY: ; Written by: Markus Fraenz,MPAE jan 94. ; 14/12/94 2.0 MeV offset for B in the upper half ; generally desactivated MF ; 18/06/95 allows for dE > A2 threshold MF,MKR ;- FUNCTION CHANTOEN,chan,xybit,hmodbit ; channel-energy-conversion: E[MeV]=(channel+facB)/facA ; facA=[128.6,2.558,24.39,12.] ; ADC1 : Erange0,Erange1,dEimod,dEhmod facB=[8.005,0.7673,7.317,7.] ; ADC1 : Erange0,Erange1,dEimod,dEhmod ; IF xybit EQ 1 THEN en=(chan+1.0+facB(2+hmodbit))/facA(2+hmodbit) >0.43 $ ELSE BEGIN nel=n_elements(chan) IF nel EQ 1 THEN BEGIN IF chan LT 256 THEN $ en=(chan+1.0+facB(0))/facA(0) $ ELSE en=(chan-256+1+facB(1))/facA(1) ;+2.0 ; with doubtfull 2.0 offset ; upper energy range ought to start at 0 again ENDIF ELSE BEGIN er0idx=where(chan LT 256,er0count) er1idx=where(chan GE 256,er1count) en=FLTARR(nel) IF er0count GT 0 THEN en(er0idx)=(chan(er0idx)+1.0+facB(0))/facA(0) IF er1count GT 0 THEN en(er1idx)= $ (chan(er1idx)-256.0+1.0+facB(1))/facA(1);+2.0 ENDELSE ENDELSE return,en end Python Code # Import Library import numpy as np import os # Declare Path TAB_path = "" # Extract the row from TAB Files def read_tab_file(filename): #Store data in list data = [] with open(filename, 'r') as file: for line in file: columns = line.strip().split() # Convert data to floats data.append([float(column) for column in columns]) # Build into a Data Matrix data = np.array(data) return data # FUNCTION CHAN -> ENERGY def chantoen(chan, xybit, hmodbit): # Formula: # E[MeV]=(channel+facB)/facA facA = np.array([128.6, 2.558, 24.39, 12.]) # ADC1: Erange0, Erange1, dEimod, dEhmod facB = np.array([8.005, 0.7673, 7.317, 7.]) # ADC1: Erange0, Erange1, dEimod, dEhmod # hmodbit must be an int if np.size(hmodbit) > 1: hmodbit = [int(x) for x in hmodbit] else: hmodbit = int(hmodbit) # Check if all elements of xybit are equal to 1 if (xybit == 1).all(): en = (chan + 1.0 + facB[2 + hmodbit]) / facA[2 + hmodbit] if en <= 0.43: en = 0.0 else: nel = np.size(chan) if nel == 1: if chan < 256: en = (chan + 1.0 + facB[0]) / facA[0] #upper energy range ought to start at 0 again else: en = (chan - 256 + 1 + facB[1]) / facA[1] + 2.0 # with doubtful 2.0 offset else: er0idx = np.where(chan < 256)[0] er1idx = np.where(chan >= 256)[0] # en=FLTARR(nel) -> creates a single-precision floating-point array en = np.zeros(nel) if np.size(er0idx) > 0: en[er0idx] = (chan[er0idx] + 1.0 + facB[0]) / facA[0] if np.size(er1idx) > 0: en[er1idx] = (chan[er1idx] - 256.0 + 1.0 + facB[1]) / facA[1] + 2.0 return en # Conversion for file in os.listdir(TAB_path): if file.endswith(".TAB"): filepath = os.path.join(TAB_path, file) data = read_tab_file(filepath) energy = chantoen(data[:,0], data[:,1], data[:,2]) print(f"{file}: {len(energy)} eV") -------------------------------------------------------------------------- EPAC_EV -------------------------------------------------------------------------- ev=fltarr(512,256) b=bytarr(3) !x.style=1 !y.style=1 loadct,40 set_plot,'ps' device,filename='epac-ev.ps',/color,xoffset=3.5,yoffset=27.0,$ xsize=26,ysize=16,/landscape filename='EV .BIN' for iyear=92,92 do begin strput,filename,string(long(iyear),format="(i2.2)"),2 for doy=25,42 do begin strput,filename,string(long(doy),format="(i3.3)"),4 openr,1,filename,error=ierror,/binary,/NONE,/BLOCK if ierror eq 0 then begin while not eof(1) do begin readu,1,b idx1=b(0) idx2=b(1) flag=b(2) idx1=idx1+255*FIX(flag AND 128b)/128 TID=FIX(flag AND 12b)/4 CeID=FIX(flag AND 3b) ev(idx1,idx2)=ev(idx1,idx2)+1 endwhile close,1 endif endfor endfor if max(ev) gt 0 then maxev=alog10(max(ev)) D1=fltarr(512) D2=fltarr(256) plot,D2,D1,xrange=[0.0,512.0],yrange=[0.0,256.0],/nodata,$ xtitle='B channel',ytitle='A channel',$ title='EPAC-EV: ' DE1=0.46 DE2=0.46 for i=0,511 do begin for j=0,255 do begin E1=float(i+0.5) E2=float(j+0.5) EE1=[E1-DE1,E1+DE1,E1+DE1,E1-DE1] EE2=[E2-DE2,E2-DE2,E2+DE2,E2+DE2] if ev(i,j) ne 0 then begin icol=fix(alog10(ev(i,j))*254/maxev) polyfill,EE1,EE2,color=icol endif endfor endfor device,/close end # Import Libraries import numpy as np import matplotlib.pyplot as plt # Initialize arrays ev = np.zeros((512, 256), dtype=float) # Process Data filename_prefix = 'EV' for iyear in range(92, 93): for doy in range(25, 43): filename = f"{filename_prefix}{iyear:02d}{doy:03d}.BIN" try: with open(filename, 'rb') as file: while True: #Read binary data byte_data = file.read(3) if not byte_data or len(byte_data) < 3: break idx1 = byte_data[0] + 255 * (byte_data[2] & 0x80) // 128 idx2 = byte_data[1] ev[idx1, idx2] += 1 except FileNotFoundError: pass # Plotting plt.figure(figsize=(26, 16)) plt.gca().set_aspect('equal', adjustable='box') plt.gca().set_xlim(0, 512) plt.gca().set_ylim(0, 256) # Get max value for scaling the visualization if np.max(ev) > 0: maxev = np.log10(np.max(ev)) else: maxev = 0 DE1 = 0.46 DE2 = 0.46 for i in range(512): for j in range(256): E1 = i + 0.5 E2 = j + 0.5 EE1 = [E1 - DE1, E1 + DE1, E1 + DE1, E1 - DE1] EE2 = [E2 - DE2, E2 - DE2, E2 + DE2, E2 + DE2] if ev[i, j] != 0: icol = int(np.fix(np.log10(ev[i, j]) * 254 / maxev)) plt.fill(EE1, EE2, color=plt.cm.jet(icol)) plt.xlabel('B channel') plt.ylabel('A channel') plt.title('EPAC-EV') plt.savefig('epac_ev.png', dpi=300, bbox_inches='tight')