#include "filter.h" #include "windmisc.h" #include "p3d_dcm.h" #include "pl_dcm.h" #include "emom_dcm.h" #include "pmom_dcm.h" #include "hkp_dcm.h" #include "eAtoD_dcm.h" #include "sst_dcm.h" #include "p3d_time.h" #include "fpc_dcm.h" #include #include #include #ifdef __cplusplus extern "C" { #endif #if 0 double str_to_time_idl(int argc,void *argv[]) { double t; IDL_STRING str; str = *(IDL_STRING *)argv[0]; if(str.length) t = YMDHMS_to_time(str.s); else t = 0; return(t); } int time_to_str_idl(int argc,void *argv[]) { double t; char *s; char *buff; /* must be at least 40 char long */ t = *(double *)argv[0]; buff = (char *)argv[1]; s = time_to_YMDHMS(t); strcpy(buff,s); return(strlen(s)); } #endif int load_data_files_idl(int argc,void *argv[]) { double *begintime,*endtime; char *mstfilename="mastfile"; /* default value */ IDL_STRING str; int ok; long size; static long memsize; static char *bigmem; if(argc<1){ printf("Must pass time interval\n"); return(0); } begintime = (double *)argv[0]; endtime = begintime + 1; if(argc>=2){ str = *((IDL_STRING *)argv[1]); if(str.length) mstfilename = str.s; } if(argc >=3){ size = *((int2 *)argv[2]); } else size=25; if(argc >= 4){ pkquality = *((int2 *)argv[3]); } else pkquality = 0; if(pkquality > 0) printf("Packet quality flag set to %d\n",pkquality); size = size *1024 * 1024; /* printf("%s (%f)\n",time_to_YMDHMS(*begintime),*begintime); */ /* printf("%s (%f)\n",time_to_YMDHMS(*endtime),*endtime); */ if(size != memsize){ if(bigmem) free(bigmem); bigmem = (char *) malloc(size); memsize = size; } if(bigmem == 0){ fprintf(debug,"Unable to allocate memory\r\n"); return(0); } ok=load_all_data_files_p(mstfilename,begintime,endtime,memsize,bigmem); if (ok > 0) { realloc((void *)bigmem,memsize - ok); memsize = memsize - ok; } return(ok); } int change_pkquality(int argc,void *argv[]) { if(argc<1){ printf("Must provide packet quality\r\n"); return(0); } pkquality = *((int2 *)argv[0]); printf("Changed packet quality to %d\r\n",pkquality); return(1); } typedef struct { int4 *options; /* options[0] is the data type, options[1] is index */ double *time; /* time[0]=sample time; time[1]=integration time */ int2 *sizes; /* 0: nphi; 1:ntheta; 2:nenergies; 3:nbins */ int2 *ptmap; /* bin map [ntheta][nphi] */ float *pt_limits; /* 0: phimin, 1:phimax, 2:thetamin; 3:thetamax */ float *data; /* counts data [nbins][nenergies] */ float *nrgs; /* energy steps [nbins][nenergies] */ float *geom; /* geometric factors [nbins] */ float *thetas; /* theta values [nbins][nenergies] */ float *phis; /* phi values [nbins][nenergies] */ float *dnrgs; /* delta energy values [nbins][nenergies] */ float *dtheta; /* delta theta values [nbins] */ float *dphi; /* delta phi values [nbins] */ float *domega; /* solid angle [nbins] */ float *eff; /* duty cycle of each energy step [nenergies] */ float *feff; /* foil electron efficiencies [nenergies] */ short *advance; /* get next point flag */ int4 *spin; /* spin */ short *magel; /* mag elevation */ short *magaz; /* mag azimuth */ } idl_3d_data; int fill_esa_data(idl_3d_data idl); int make_ptmap(data_map_3d *map,idl_3d_data idl); int fill_pl_data(idl_3d_data idl); int fill_plb_data(idl_3d_data idl); int fill_sst_foil_data(idl_3d_data idl); int fill_sst_open_data(idl_3d_data idl); int fill_sst_ft_data(idl_3d_data idl); int fill_sst_ot_data(idl_3d_data idl); int fill_sst_f_spectra_data(idl_3d_data idl); int fill_sst_o_spectra_data(idl_3d_data idl); int fill_sst_t_spectra_data(idl_3d_data idl); int fill_sst_fast_rates(idl_3d_data idl); int fill_sst_slow_rates(idl_3d_data idl); /* returns 0 on error */ int get_3dbins_idl(int argc,void *argv[]) { idl_3d_data idl; /* contains all input/output pointers */ int data_type; /*printf("argc= %d\n",argc); */ idl.options = (int4 *)argv[0]; /*options[0] is data type; *options[1] is index, -1 == get by time *options[2] is selection mode for sst*/ idl.time = (double *)argv[1]; idl.sizes = (int2*)argv[2]; idl.ptmap = (int2 *)argv[3]; idl.pt_limits = (float *)argv[4]; idl.data = (float *)argv[5]; idl.nrgs = (float *)argv[6]; idl.geom = (float *)argv[7]; idl.thetas = (float *)argv[8]; idl.phis = (float *)argv[9]; idl.dnrgs = (float *)argv[10]; idl.dtheta = (float *)argv[11]; idl.dphi = (float *)argv[12]; idl.domega = (float *)argv[13]; idl.eff = (float *)argv[14]; idl.feff = (float *)argv[15]; idl.spin = (int4 *)argv[16]; idl.magel = (int2 *)argv[17]; idl.magaz = (int2 *)argv[18]; idl.advance = (int2 *)argv[19]; data_type = idl.options[0]; switch(data_type){ case 0: /* EH */ case 1: /* EL */ case 2: /* PH */ return( fill_esa_data(idl) ); case 3: /* PL */ return( fill_pl_data(idl) ); case 4: /* SST FOIL */ case 19: /* SST FOIL BURST*/ return( fill_sst_foil_data(idl) ); case 5: /* SST OPEN */ case 20: /* SST OPEN BURST*/ return( fill_sst_open_data(idl) ); case 7: /* PL Burst */ return( fill_plb_data(idl) ); case 8: /* SST F+T */ return( fill_sst_ft_data(idl) ); case 9: /* SST O+T */ return( fill_sst_ot_data(idl) ); case 10: /* SST FOIL SPECTRA */ case 17: /* SST FOIL BURST SPECTRA */ return( fill_sst_f_spectra_data(idl) ); case 11: /* SST OPEN SPECTRA */ case 18: /* SST OPEN BURST SPECTRA */ return( fill_sst_o_spectra_data(idl) ); case 12: /* SST F+T % O+T SPECTRA */ return( fill_sst_t_spectra_data(idl) ); case 13: /* EL Burst */ return( fill_esa_data(idl) ); case 14: /* EH Slice */ return( fill_esa_data(idl) ); case 15: /* PH Burst */ return( fill_esa_data(idl) ); case 16: /* Eesa Low Cuts */ return( fill_esa_data(idl) ); case 21: /* SST FAST RATES */ return( fill_sst_fast_rates(idl) ); case 22: /* SST SLOW RATES */ return( fill_sst_slow_rates(idl) ); default: return(0); } } int get_fpc_idl(int argc,void *argv[]) { /* argv[0] will be the fpc_data struct */ fpc_xcorr_str *fpc_idl = (fpc_xcorr_str *) argv[0]; int advance = *((short *) argv[1]); static packet_selector pks; if (advance) { SET_PKS_BY_INDEX(pks,pks.index+advance,FPC_D_ID); } else if (fpc_idl->select_by == BY_INDEX) { SET_PKS_BY_INDEX(pks,fpc_idl->index,FPC_D_ID); } else { SET_PKS_BY_TIME(pks,fpc_idl->time,FPC_D_ID); } /* fill idl struct */ return get_next_fpc(&pks, fpc_idl); } int f_bin_channels[48] = {3,3,2,2,2,2,4,4,4,4,4,4,4,4, 5,5,5,5,0,0,1,1,1,1,3,3,2,2, 2,2,4,4,4,4,4,4,4,4,5,5,5,5, 0,0,1,1,1,1}; int o_bin_channels[48] = {0,0,5,5,5,5,4,4,4,4,4,4,4,4,2,2, 2,2,3,3,1,1,1,1,0,0,5,5,5,5,4,4, 4,4,4,4,4,4,2,2,2,2,3,3,1,1,1,1}; int fill_sst_foil_data(idl_3d_data idl) { static sst_3d_F_distribution sst; double time; int e,p,t,bin,c; int ne,np,nt,nb; int ok,cntr; uint2 validmask; packet_selector *pksp; static packet_selector pks_sf; static packet_selector pks_sfb; int id; switch( idl.options[0] ){ case 4: /* SF */ pksp = &pks_sf; id = S_3D_F_ID; break; case 19: /* SFB */ pksp = &pks_sfb; id = S_HS_BST_ID; break; default: ok = 0; } if ( *idl.advance ) { SET_PKS_BY_INDEX(*pksp,pksp->index+*idl.advance,id) ; } else if (idl.options [1] < 0) { /* negitive option[1] means get by time*/ SET_PKS_BY_TIME(*pksp,idl.time[0],id) ; } else { SET_PKS_BY_INDEX(*pksp,idl.options[1],id) ; } ok = get_next_sst_3d_F_str(pksp,&sst, idl.options[2], &validmask); if(!ok) return(0); idl.options[1] = pksp->index; idl.time[0] = sst.time; idl.time[1] = sst.integ_t; idl.time[2] = sst.mass; idl.time[3] = sst.geom_factor; np = idl.sizes[0] = 32; nt = idl.sizes[1] = 15; ne = idl.sizes[2] = 7; nb = idl.sizes[3] = 48; idl.pt_limits[0] = -90.; idl.pt_limits[1] = -180.; idl.pt_limits[2] = 90.; idl.pt_limits[3] = 180.; cntr = 0; for(t=0;t>e)) ; break; case 19: /* SFB */ for(e=0;eindex+*idl.advance,id) ; } else if (idl.options [1] < 0) { /* negitive option[1] means get by time*/ SET_PKS_BY_TIME(*pksp,idl.time[0],id) ; } else { SET_PKS_BY_INDEX(*pksp,idl.options[1],id) ; } ok = get_next_sst_3d_O_str(pksp,&sst, idl.options[2], &validmask); if(!ok) return(0); idl.options[1] = pksp->index; idl.time[0] = sst.time; idl.time[1] = sst.integ_t; idl.time[2] = sst.mass; idl.time[3] = sst.geom_factor; np = idl.sizes[0] = 32; nt = idl.sizes[1] = 15; ne = idl.sizes[2] = 9; nb = idl.sizes[3] = 48; idl.pt_limits[0] = -90.; idl.pt_limits[1] = -180.; idl.pt_limits[2] = 90.; idl.pt_limits[3] = 180.; cntr = 0; for(t=0;t>e)); break; case 20: /* SOB */ for(e=0;eindex+*idl.advance,id) ; } else if (idl.options [1] < 0) { /* negitive option[1] means get by time*/ SET_PKS_BY_TIME(*pksp,idl.time[0],id) ; } else { SET_PKS_BY_INDEX(*pksp,idl.options[1],id) ; } ok = get_next_sst_3d_O_str(pksp,&sst, idl.options[2], &validmask); if(!ok) return(0); idl.options[1] = pksp->index; idl.time[0] = sst.time; idl.time[1] = sst.integ_t; idl.time[2] = sst.mass; idl.time[3] = sst.geom_factor; np = idl.sizes[0] = 32; nt = idl.sizes[1] = 15; ne = idl.sizes[2] = 1; nb = idl.sizes[3] = 14; idl.pt_limits[0] = -90.; idl.pt_limits[1] = -180.; idl.pt_limits[2] = 90.; idl.pt_limits[3] = 180.; cntr = 0; for(t=0;tindex+*idl.advance,id) ; } else if (idl.options [1] < 0) { /* negitive option[1] means get by time*/ SET_PKS_BY_TIME(*pksp,idl.time[0],id) ; } else { SET_PKS_BY_INDEX(*pksp,idl.options[1],id) ; } ok = get_next_sst_spectra_str(pksp,&sst); if(!ok) return(0); idl.options[1] = pksp->index; idl.time[0] = sst.time; idl.time[1] = sst.integ_t; idl.time[2] = sst.mass; idl.time[3] = sst.geom_factor; np = idl.sizes[0] = 32; nt = idl.sizes[1] = 15; ne = idl.sizes[2] = 16; nb = idl.sizes[3] = 6; idl.pt_limits[0] = -90.; idl.pt_limits[1] = -180.; idl.pt_limits[2] = 90.; idl.pt_limits[3] = 180.; cntr = 0; for(t=0;tindex+*idl.advance,id) ; } else if (idl.options [1] < 0) { /* negitive option[1] means get by time*/ SET_PKS_BY_TIME(*pksp,idl.time[0],id) ; } else { SET_PKS_BY_INDEX(*pksp,idl.options[1],id) ; } ok = get_next_sst_spectra_str(pksp,&sst); if(!ok) return(0); idl.options[1] = pksp->index; idl.time[0] = sst.time; idl.time[1] = sst.integ_t; idl.time[2] = sst.mass; idl.time[3] = sst.geom_factor; np = idl.sizes[0] = 32; nt = idl.sizes[1] = 15; ne = idl.sizes[2] = 24; nb = idl.sizes[3] = 6; idl.pt_limits[0] = -90.; idl.pt_limits[1] = -180.; idl.pt_limits[2] = 90.; idl.pt_limits[3] = 180.; cntr = 0; for(t=0;t1) { idl.nrgs[cntr] = sst.e_OT_mid[bin-2][e]; idl.dnrgs[cntr] = sst.e_OT_max[bin-2][e]-sst.e_OT_min[bin-2][e];} cntr++; } idl.dtheta[bin] = sst.dtheta[bin]; idl.dphi[bin] = sst.dphi[bin]; idl.geom[bin] = sst.geom[bin]; idl.domega[bin] = sst.domega[bin]; } for(e=0;eindex+*idl.advance,id) ; } else if (idl.options [1] < 0) { /* negitive option[1] means get by time*/ SET_PKS_BY_TIME(*pksp,idl.time[0],id) ; } else { SET_PKS_BY_INDEX(*pksp,idl.options[1],id) ; } ok = get_next_p3d(pksp, &map, exppk); make_ptmap(&map,idl); idl.time[0] = map.time; idl.time[1] = map.integ_t; idl.time[2] = map.mass; idl.time[3] = map.geom_factor; *idl.spin = map.spin; idl.options[1] = pksp->index; return(ok); } static int tsect[32] = {0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3, 4, 5, 6, 7, 8, 9,10,11,12,12,13,13,14,14,14,14,15,15,15,15}; int make_ptmap(data_map_3d *map,idl_3d_data idl) { int p,t,e; int np,nt,ne,nb; int ts,n,offset,b; int shift,ps; int cntr; double flux; /* Map stuff: */ np = idl.sizes[0] = 32; nt = idl.sizes[1] = 32; ne = idl.sizes[2] = map->nenergies; nb = idl.sizes[3] = map->nbins; shift = map->shift; for(t=0;t<32;t++){ ts = tsect[t]; for(p=0;p<32;p++){ b = map->ptmap[ts*32+p]; ps = (p - map->p0) & 31; idl.ptmap[t*np+ps] = b; /* #define MAKE_PTMAP_DEBUG*/ #if defined (MAKE_PTMAP_DEBUG) printf("t: %d, p: %d, ts: %d, ps: %d, idl.ptmap[%d]: %d\r\n", t,p,ts,ps,t*np+ps,idl.ptmap[t*np+ps]); #endif defined (MAKE_PTMAP_DEBUG) } } idl.pt_limits[0] = -90.; idl.pt_limits[1] = -180.; idl.pt_limits[2] = 90.; idl.pt_limits[3] = 180.; /* data and energy steps: */ *idl.magel = map->magel; *idl.magaz = map->magaz; cntr =0; for(b=0;bnbins;b++){ idl.geom[b] = map->bin[b].geom; idl.dtheta[b] = map->bin[b].dtheta; idl.dphi[b] = map->bin[b].dphi; idl.domega[b] = map->bin[b].domega; offset = map->bin[b].offset; n = map->bin[b].ne; for(e=0;edata[offset+e]; } else{ /* twice as many e-steps as overall map e-steps */ flux = map->data[offset+2*e]; flux += map->data[offset+2*e+1]; } idl.data[cntr] = flux; if(ne == 30){ idl.nrgs[cntr] = map->nrg30[e]; idl.dnrgs[cntr] = map->dnrg30[e]; } else { idl.nrgs[cntr] = map->nrg15[e]; idl.dnrgs[cntr] = map->dnrg15[e]; } idl.thetas[cntr] = map->bin[b].theta; idl.phis[cntr] = map->bin[b].phi[e]; cntr++; } } for(e=0;es, "el") ) pkid = E3D_88_ID ; else if ( !strcmp (type->s, "eh") ) pkid = E3D_UNK_ID ; else if ( !strcmp (type->s, "pl") ) pkid = PLSNAP_ID ; else if ( !strcmp (type->s, "ph") ) pkid = P3D_ID ; else if ( !strcmp (type->s, "fpc") ) pkid = FPC_D_ID ; else if ( !strcmp (type->s, "sf") ) pkid = S_3D_F_ID ; else if ( !strcmp (type->s, "so") ) pkid = S_3D_O_ID ; else if ( !strcmp (type->s, "fr") ) pkid = S_3D_O_ID ; else if ( !strcmp (type->s, "sft") ) pkid = S_T_DST_ID ; else if ( !strcmp (type->s, "sot") ) pkid = S_T_DST_ID ; else if ( !strcmp (type->s, "fspc") ) pkid = S_RATE3_ID ; else if ( !strcmp (type->s, "ospc") ) pkid = S_RATE3_ID ; else if ( !strcmp (type->s, "tspc") ) pkid = S_RATE3_ID ; else if ( !strcmp (type->s, "sr") ) pkid = S_RATE3_ID ; else if ( !strcmp (type->s, "plb") ) pkid = P_SNAP_BST_ID ; else if ( !strcmp (type->s, "elb") ) pkid = E3D_BRST_ID; else if ( !strcmp (type->s, "ehs") ) pkid = FPC_P_ID; else if ( !strcmp (type->s, "phb") ) pkid = P3D_BRST_ID; else if ( !strcmp (type->s, "elc") ) pkid = E3D_CUT_ID; else if ( !strcmp (type->s, "sob") ) pkid = S_HS_BST_ID; else if ( !strcmp (type->s, "sfb") ) pkid = S_HS_BST_ID; else { fprintf(stderr, "get_time_array: unknown type\r\n"); return 0; } /* get the time array now */ return get_time_points(pkid, *max_array, time_array) ; } #ifdef __cplusplus } /* End extern C */ #endif