/* * Combined STARE analysis/display program PJ Oct-23 97 Last mod Nov-12 98 * --------------------------------------- * * The program operates in three modes. * * (1) Raw data -> analyzed data converter mode * This mode is selected if Tela was called with the -b * flag (BatchMode() is nonzero) AND variable 'inputfile' has been set * a string value (using Tela's -e option, see tela --help). * Typical usage: tela -be 'inputfile="02072000.n97"' stare.t * This produces file 02072000.n97.analyzed.hdf.gz. It is old-style * HDF file with 32-bit floats. Quantization for snr/vel is used so * that gzip can compress it down to about 250KB. * If you don't mind any messages, add the -s option for Tela. * * (2) Analyzed data -> RTI plot mode * This mode is selected if Tela is in batch mode AND * variable 'analyzedfile' has been set to a string value. * The analyzedfile must be a full path name relative to /proj/STARE/analyzed, * for example: * tela -be 'analyzedfile="1997/12/19971202.fin.analyzed.hdf.gz"' stare.t * This command writes /proj/STARE/RTI/1997/12/19971202.RTI.gif. * * (3) Analyzed data -> LTV plot mode. * This mode is selected if Tela is in batch mode AND * variable 'ltvfile' has been set to a string value. * The parameter duration, integration and ltv_dir can also be set beforehand: * - 'duration' is the length of the event in minutes (default 1440, i.e. 24 hours) * - 'integration' is the integration step in seconds (default 600, i.e. 10 minutes) * - ltv_dir is a string * Otherwise similar as mode (2). * * (4) Interactive mode: * * Just type tela stare.t * * It looks for analyzed files under /proj/STARE/analyzed and tries * to load the latest day file. Other files can be loaded by "Load new day file ...". * It looks for both uncompressed and gzipped analyzed files. * Gzipped files are uncompressed in $TMPDIR, so write permission to /proj/STARE/analyzed * is not needed. The file is removed from $TMPDIR immediately after loading it in. */ longcoeff = 0.45; lat1 = 67; lat2 = 72; long1 = 14; long2 = 27; lat1_map = 60; lat2_map = 73; long1_map = 4; long2_map = 30; lat_hankasalmi = 62.3047; long_hankasalmi = 26.6494; lat_malvik = 63.6667; long_malvik = 10.73; doplot = 1; // always 1, except temporarily in LTVplot angles_malvik = (pi/180)*(26.2+(-3.5:3.5)*3.6); angles_hankasalmi = (pi/180)*(-15.5 + (-3.5:3.5)*3.6); NorwegianSign = -1; FinnishSign = +1; monthstr = strmat("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"); time_wanted = 0; DisplayMap = on; DisplayBeams = on; DisplayFinnish = on; DisplayNorwegian = on; DisplayVectors = on; DisplayBothVectors = on; DisplaySNR = on; DisplayRTI = off; PlotEField = off; UseInterpolation = off; UseColorPrinting = on; DBLIMIT = 5.0; UseStrictDBLimit = on; // If UseStrictDBLimit is nonzero, both radars must have their SNR above DBLIMIT for a vector to be shown. // If UseStrictDBLimit is zero, it suffices if the average of F and N is above DBLIMIT. DataDir = "/proj/STARE/analyzed"; if (BatchMode() && isstring(inputfile)) { BatchAnalysisOperation(inputfile); return; } else if (BatchMode() && isstring(analyzedfile)) { RTIBatchAnalysisOperation(analyzedfile); return; } else if (BatchMode() && isstring(ltvfile)){ if (!isreal(duration)) { duration = 1440; format("No duration supplied. Setting it to default value ``\n",duration); }; if (!isstring(ltv_dir)) { ltv_dir = sformat("/tmp/``/miracle",getenv("USER")); format("No ltv tmp-directory (ltv_dir) supplied. Setting it to default value ``\n",ltv_dir); }; if (!isreal(integration)) { integration = 600; format("No integration supplied. Setting it to default value ``\n",integration); }; format("LTV batch mode: duration is ``\n",duration); format("LTV batch mode: integration is ``\n",integration); LTVBatchAnalysisOperation(ltvfile,duration,integration,ltv_dir); return; }; UpdateColorMode(); //cd(DataDir); currentyear = run("date +%Y"); currentyear = currentyear[1:4]; M = strmat2(run(sformat("cd ``; find `` -name '*.analyzed.hdf.gz' -print | sort",DataDir,currentyear))); [Nfiles,dum] = size(M); if (Nfiles > 0) { year = evalexpr(M[Nfiles,1:4]); month = evalexpr(M[Nfiles,6:7]); day = evalexpr(M[Nfiles,15:16]); format("Starting with latest file found under ``\n",DataDir); } else { format("*** No files found in current directory, trying 19971010 ...\n"); year = 1997; month = 10; day = 9 + 1; }; success = ReadDataFiles(); while (!success) { format("Could not read last data file (``````).\n",year,month,day); [year,month,day,cancelled] = AskDay(); if (cancelled) quit(); }; figure(1); UpdatePlot(); stare(); function stare() global { repeat ans = smenu("Choose function:", "Load new day file ...", "Next time", "Previous time", "One minute forward", "One minute backward", "5 minutes forward", "5 minutes backward", "One hour forward", "One hour backward", "Goto given time ...", "Toggles submenu...", "Make animation ...", "Make LTV plot and export to file ...", "Save previous snr&vectors to HDF file", "Save previous snr&vectors to ASCII file", "Print data values on given lat,long point", "Quit"); if (strstarteq(ans,"Load new day")) { [year,month,day,cancelled] = AskDay(); if (!cancelled) { success = ReadDataFiles(); time_wanted = 0; UpdatePlot(); if (DisplayRTI) { figure(2); RTIplots(); figure(1); }; }; } else if (strstarteq(ans,"Next")) { time_wanted = limit(time_wanted+integration_time,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"Previous")) { time_wanted = limit(time_wanted-integration_time,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"One minute forward")) { time_wanted = limit(time_wanted+60,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"One minute backward")) { time_wanted = limit(time_wanted-60,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"5 minutes forward")) { time_wanted = limit(time_wanted+5*60,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"5 minutes backward")) { time_wanted = limit(time_wanted-5*60,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"One hour forward")) { time_wanted = limit(time_wanted+3600,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"One hour backward")) { time_wanted = limit(time_wanted-3600,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"Goto given time")) { hour = input("Hour"); minute = input("Minute"); time_wanted = limit(3600*hour+60*minute,0,24*3600); UpdatePlot(); } else if (strstarteq(ans,"Toggle")) { ToggleMenu(); } else if (strstarteq(ans,"Make anim")) { MakeAnimation(); } else if (strstarteq(ans,"Make LTV")) { figure(3); LTVplot(); figure(1); } else if (strstarteq(ans,"Save previous snr&vectors to HDF")) { GrabLastToHdfFile(); } else if (strstarteq(ans,"Save previous snr&vectors to ASCII")) { GrabLastToAsciiFile(); } else if (strstarteq(ans,"Print data values")) { atlat = input("Latitude "); atlon = input("Longitude"); PrintDataAt(atlat,atlon); } else if (strstarteq(ans,"Quit")) { break; } else { format("*** unknown menu selection, internal error?\n"); }; until 0; }; // end of function stare() function ToggleMenu() local(ans) { repeat ans = smenu("Choose:", sformat("Show RTI plots (now ``)",where(DisplayRTI,"on","off")), sformat("Toggle map display (now ``)",where(DisplayMap,"on","off")), sformat("Toggle beam display (now ``)",where(DisplayBeams,"on","off")), where(successF, sformat("Toggle Finnish data display (now ``)",where(DisplayFinnish,"on","off")), " (Toggle Finnish data display, off since no data)"), where(successN, sformat("Toggle Norwegian data display (now ``)",where(DisplayNorwegian,"on","off")), " (Toggle Norwegian data display, off since no data)"), sformat("Toggle vector display (now ``)",where(DisplayVectors,"on","off")), sformat("Toggle display of both components always (now ``)",where(DisplayBothVectors,"on","off")), sformat(where(PlotEField,"Show electron flow (now electric field)","Show electric field (now electron flow)")), sformat("Toggle SNR display (now ``)",where(DisplaySNR,"on","off")), sformat("Toggle 'strict' way of dB limiting (now ``)",where(UseStrictDBLimit,"on","off")), sformat("Set threshold dB value for vector plotting (now `` dB)",DBLIMIT), sformat("Toggle pcolor interpolation (now ``)",where(UseInterpolation,"on","off")), sformat("Toggle color PostScript (now ``)",where(UseColorPrinting,"on","off")), "Quit this menu"); if (strstarteq(ans,"Show RTI")) { DisplayRTI = !DisplayRTI; if (DisplayRTI) { figure(2); RTIplots(); figure(1); } else { closefig(2); }; } else if (strstarteq(ans,"Show ele")) { PlotEField = !PlotEField; UpdatePlot(); } else if (strstarteq(ans,"Toggle map")) { DisplayMap = !DisplayMap; UpdatePlot(); } else if (strstarteq(ans,"Toggle beam")) { DisplayBeams = !DisplayBeams; UpdatePlot(); } else if (strstarteq(ans,"Toggle Finnish")) { DisplayFinnish = !DisplayFinnish; UpdatePlot(); } else if (strstarteq(ans,"Toggle Norwegian")) { DisplayNorwegian = !DisplayNorwegian; UpdatePlot(); } else if (strstarteq(ans,"Toggle vector")) { DisplayVectors = !DisplayVectors; UpdatePlot(); } else if (strstarteq(ans,"Toggle display of both")) { DisplayBothVectors = !DisplayBothVectors; UpdatePlot(); } else if (strstarteq(ans,"Toggle SNR")) { DisplaySNR = !DisplaySNR; UpdatePlot(); } else if (strstarteq(ans,"Set threshold")) { if (UseStrictDBLimit) format("The criterion is that both radar SNRs must be larger than DBLIMIT.\n") else format("The criterion is that the arithmetic average of the radar SNRs must be larger than DBLIMIT.\n"); format("Currently DBLIMIT is set to `` dB\n",DBLIMIT); DBLIMIT = input("Give new DBLIMIT "); UpdatePlot(); } else if (strstarteq(ans,"Toggle 'strict'")) { UseStrictDBLimit = !UseStrictDBLimit; if (UseStrictDBLimit) format("Both radars must have SNR > `` dB for a vector to be shown\n") else format("Average SNR of N and F radars must be > `` dB for a vector to be shown\n"); UpdatePlot(); } else if (strstarteq(ans,"Toggle pcolor interpolation")) { UseInterpolation = !UseInterpolation; UpdatePlot(); } else if (strstarteq(ans,"Toggle color")) { UseColorPrinting = !UseColorPrinting; UpdateColorMode(); format("Closing and reopening figure(s) ...\n"); closefig(1); // force PlotMTV to restart UpdatePlot(); if (DisplayRTI) { closefig(2); figure(2); RTIplots(); figure(1); }; } else if (strstarteq(ans,"Quit")) { break; } else { format("*** unknown menu selection, internal error?\n"); }; until 0; }; // end of function ToggleMenu // ------------- Some Matlab style functions ----------- function p = angle(h) {p = Im(log(h))}; // ----------------------------------------------------- function [year,month,day,cancelled] = AskDay() { current_year = evalexpr(run("date +%Y")); s = input_string(sformat("Year [``, c to cancel]",current_year)); if (!streq(s,"c")) { if (length(s)==0) year = current_year else year = evalexpr(s); month = input("Month (1-12)"); day = input("Day of month (1-31)"); cancelled = 0; } else cancelled = 1; }; function y = readlong(fp) { y = 0; i1 = fgetc(fp); i2 = fgetc(fp); i3 = fgetc(fp); i4 = fgetc(fp); if (feof(fp)) return; y = 16777216*i4 + 65536*i3 + 256*i2 + i1; }; /* * Data file format * ---------------- * The data file consists of 20-second data dumps (which we shall call records). * Each record contains len 32-bit integers data[i]. * data[1] is equal to (len-1)*4, i.e., the number of bytes in data[2:len]. * data[2:101] is the parameter block (parbl). * data[102:len] is the measured data itself. * Each record is independent of the other records. */ function [parbl,data,success] = read_raw_1(fid,len) // Read the next 20-second dump (record) from previously opened file // with file descriptor fid. The length of the record must be known and // passed in len (the number of 32-bit integers in the record). // success is set to 0 on error, 1 on normal completion, and 2 on EOF. { success = 1; cdata = fread(fid,4*len); // Read len 32-bit integers CL = length(cdata); if (CL != 4*len) { // Give an error message if some data could be read, but not all. // If no data could be read, we are most probably at the end of the file, // and no error message needs to be given. if (CL > 0) { format("*** stare.t:read_raw_1(): Read `` bytes instead of ``\n", CL,4*len); success = 0; } else success = 2; return; }; // Convert character data cdata into 32-bit integer data data = cdata[1:4:CL]; data+= 256*cdata[2:4:CL]; data+= 65536*cdata[3:4:CL]; data+= 16777216*cdata[4:4:CL]; parbl = data[2:101]+0; // add 0 to get rid of possible strflag data = data[102:len]+0; // also }; /* function [parbl,data,success] = read_raw(inputfile) // Read the whole raw input file. Return parameter blocks in parbl[1:100,1:n] // and data itself in data[:,1:n], where n is the number of 20-second dumps (records) // in the file. All records are read by this function. // success is set to 1 if no error occurred, otherwise to 0. global(fread) { fid = fopen(inputfile,"r"); success = 1; if (fid < 0) {success = 0; return}; len = round(readlong(fid)/4) + 1; if (feof(fid)) {success = 0; return}; fclose(fid); fid = fopen(inputfile,"r"); // MAXDATA = 800000; // Maximum number of long ints (long=4 bytes) in the file MAXDATA = 19000000; // Should be enough for day-files // MAXDATA = 16800000; if (isfunction(fread)) { cdata = fread(fid,4*MAXDATA); CL = length(cdata); data = cdata[1:4:CL]; data+= 256*cdata[2:4:CL]; data+= 65536*cdata[3:4:CL]; data+= 16777216*cdata[4:4:CL]; cdata = 0; } else { format("read_raw: Using slow file read (no fread in this Tela kernel)\n"); data = izeros(MAXDATA); i = 1; while (1) { data[i] = readlong(fid); if (feof(fid)) break; i++; }; format("read `` longs\n",i-1); data = data[1:i-1]; }; N = round(length(data)/len); // format("MAXDATA=``, CL=``, N=``, len=``, length(data)=``\n",MAXDATA,CL,N,len,length(data)); data = transpose(reshape(data,N,len)); parbl = data[2:101,:]; data = data[102:len,:]; fclose(fid); }; */ function [range,sttime,entime,snr,vel,success] = analys(parbl,data) // range will be vector of length 50 // sttime,entime will be vectors of length n // snr,vel will be 3D tables (8,50,n) global(mean) { success = 1; // Store the end time of each integration in seconds // from the beginning of UT day // Store the end time of each integration in seconds // from the beginning of UT day entime = 3600.0*parbl[7,:] + 60.0*parbl[8,:] + parbl[9,:]; // Obtain the start time by subtracting the integration time sttime = entime - parbl[10,:]; Ndumps = length(entime); offset = 44; // Store the power profile control parameters to variables // the cryptic _XXXX part refers to STARE documentation PP_INTC = parbl[offset+2,1]; // Number of repetitions PP_DRAN = round(parbl[offset+5,1]/10); // Range step PP_SRAN = round(parbl[offset+6,1]/10); // First range PP_NGATS = parbl[offset+7,1]; // Number of signal gates PP_NCH = parbl[offset+8,1]; // Number of channels (beams) PP_NGATB = parbl[offset+11,1]; // Number of background gates // Sanity check if (PP_NGATS != 50 || PP_NCH != 8) { format("*** stare.t:analys: PP_NGATS=``, PP_NCH=``\n",PP_NGATS,PP_NCH); format(" File seems to be corrupt, record ignored\n"); success = 0; return; }; // Extract power profile data from the data dump ra_last = 0; ppdata = transpose(reshape(data[ra_last+(1:PP_NGATS*PP_NCH),:],PP_NGATS,PP_NCH,Ndumps),#(2,1,3)); // Last result memory address used ra_last += PP_NGATS*PP_NCH; // Extract background data from the data dump backgr = transpose(reshape(data[ra_last+(1:PP_NGATB*PP_NCH),:],PP_NGATB,PP_NCH,Ndumps),#(2,1,3)); // Average the background gates s = size(backgr); s[2] = 1; backgr = reshape(mean(backgr,2),s); ra_last += PP_NGATB*PP_NCH; // One way of calculating the signal-to-noise ratio pp = (ppdata/backgr[:,ones(PP_NGATS),:]-1); // snr snr = 10*log(max(pp,0.01))/log(10); // snr in dB // another way of doing the job (by Erling) // ppEN = (ppdata-backgr(pp_CH,:)).*10 .^(3-parbl(pp_ch+13,:))/PP_INTC; // snr = 10*log10(max(ppEN,0.01))-43; range = 0.0 + PP_SRAN + (0:PP_NGATS-1)*PP_DRAN; // Velocity offset = 56; // Store the double pulse control parameters to variables // the cryptic _XXXX part refers to STARE documentation DP_INTC = parbl[offset+2,1]; // Number of repetitions DP_LTAU = parbl[offset+3,1]; // lag resolution DP_NLAGS = parbl[offset+4,1]; // number of lags DP_DRAN = round(parbl[offset+5,1]/10); // Range step DP_SRAN = round(parbl[offset+6,1]/10); // First range DP_NGATS = parbl[offset+7,1]; // Number of signal gates DP_NCH = parbl[offset+8,1]; // Number of channels (beams) double = transpose(reshape(data[ra_last+(1:2:2*DP_NGATS*DP_NCH),:] + 1i*data[ra_last+(2:2:2*DP_NGATS*DP_NCH),:], DP_NGATS,PP_NCH,Ndumps),#(2,1,3)); ra_last += DP_NGATS*DP_NCH*2; vel_to_phase=4*pi*140e6*DP_LTAU*1e-6/3e8; vel = -angle(double)/vel_to_phase; }; function [range,sttime,entime,snr,vel,success] = read_raw_and_analys(inputfile) { fid = fopen(inputfile,"r"); success = 1; if (fid < 0) {success = 0; return}; len = round(readlong(fid)/4) + 1; if (feof(fid)) {success = 0; return}; fclose(fid); fid = fopen(inputfile,"r"); nmax = 4320; sttime = zeros(nmax); entime = sttime; snr = zeros(8,50,nmax); vel = snr; AnalysErrors = 0; for (n=1; n<=nmax; n++) { [parbl,data,success] = read_raw_1(fid,len); if (success != 1) break; // success=0: error, success=1: ok, success=2, end of file parblM = reshape(parbl,length(parbl),1); dataM = reshape(data,length(data),1); [range,sttime1,entime1,snr1,vel1,successAnalys] = analys(parblM,dataM); if (successAnalys == 1) { sttime[n] = sttime1[1]; entime[n] = entime1[1]; snr[:,:,n] = snr1[:,:,1]; vel[:,:,n] = vel1[:,:,1]; } else { AnalysErrors++; n--; if (AnalysErrors > 10) { format("*** More than 10 analysis errors, not reading file any more\n"); success = 1; break; } } }; if (success == 2) success = 1; // eof is no error n--; if (n < nmax) { format("Truncating sttime,entime,snr,vel from `` --> ``\n",nmax,n); sttime = sttime[1:n]; entime = entime[1:n]; snr = snr[:,:,1:n]; vel = vel[:,:,1:n]; if (n > 0) success = 1; // this line added 980624 PJ }; fclose(fid); }; function [x] = Quantize(coeff) { rcoeff = 1.0*coeff; // x = rcoeff*round(x/rcoeff); // Try to conserve memory as much as possible x/= rcoeff; ix = round(x); x = rcoeff*ix; }; function LTVBatchAnalysisOperation(ltvfile, duration,integration,ltv_dir) global(year,month,day,DataDir,success) { /* * tela -be 'ltvfile="1997/12/971202.fin.analyzed.hdf.gz"' stare.t * This command writes /proj/STARE/RTI/1997/12/19971202.RTI.gif. * integration, duration and the save directory can be supplied in * the option-line */ year = evalexpr(ltvfile[1:4]); month = evalexpr(ltvfile[6:7]); day = evalexpr(ltvfile[15:16]); format("year=``, month=``, day=``\n",year,month,day); cd(DataDir); success = ReadDataFiles(); if (!success) return; plotopt("-bg white -fg black -colorps -o dataplot.ps -noxplot -print -printcmd echo -Pnull"); // the following function is a rewrite of plot_LTV dump_LTV(duration,integration,ltv_dir); }; function BatchAnalysisOperation(inputfile) global(range,sttime,entime,snr,vel) { if (!SilentMode()) format("stare.t:BatchAnalysisOperation: reading/analyzing raw file \"``\"\n",inputfile); [range,sttime,entime,snr,vel,success] = read_raw_and_analys(inputfile); if (!success) {format("*** BatchAnalysisOperation: could not read input file \"``\"\n",inputfile); return}; outputfile = #(inputfile,".analyzed.hdf"); // Make 97mmdd etc. filenames -> 1997mmdd filenames if (streq(outputfile[1:2],"97") || streq(outputfile[1:2],"98") || streq(outputfile[1:2],"99")) outputfile = #("19",outputfile); if (!SilentMode()) format("stare.t:BatchAnalysisOperation: Writing file ``\n",outputfile); // Save snr and vel arrays with fewer digits. SNR is saved with 0.1 accuracy and vel 1 m/s accuracy. // This trick enables gzip to compress the file much more than otherwise: // quantization size of .analyzed.hdf.gz file // ------------ ----------------------------- // without quantize: 500 KB // quantize 0.1, 1: 240 KB // quantive 1, 10: 150 KB // Thus a further 40% reduction in file size could be achieved by storing SNR with accuracy 1 // and vel with accuracy 10 m/s, but we don't go this far now. [snr] = Quantize(0.1); [vel] = Quantize(1); DBlimit = 0.0; // A further reduction in file size is obtained if we replace values smaller than DBlimit by DBlimit. // These are essentially noise. ind = find(snr < DBlimit); vel[ind] = 0; snr[ind] = DBlimit; ind = 0; // save a little bit of memory save(outputfile,"range","sttime","entime","snr","vel"); if (!SilentMode()) format("stare.t:BatchAnalysisOperation: Compressing it --> ``.gz\n",outputfile); system(sformat("gzip -f -9 ``",outputfile)); // -f == force compression even if the file already exists // since this is batch mode function, we don't want gzip to ask any questions }; // '`´ // ------------------------------- Graphics ------------------------------- function loadmapfile() global(MAP,lat1_map,lat2_map,long1_map,long2_map, lat_malvik,long_malvik, lat_hankasalmi,long_hankasalmi) { load("MAAMERI.HDF"); ind = find(lat1_map <= MAP[:,1] && MAP[:,1] <= lat2_map && -long2_map <= MAP[:,2] && MAP[:,2] <= -long1_map || MAP[:,1] == -999 && MAP[:,2] == -999); // Remove adjacent -999 entries // format("Before remove: length(ind) = ``\n",length(ind)); n = length(ind); ind2 = ind; badvec = #(-999,-999); for ({i=2;j=2}; i<=n; i++) if (!(MAP[ind[i],:] == badvec && MAP[ind[i-1],:] == badvec)) { ind2[j] = ind[i]; j++; }; ind = ind2[1:min(j,length(ind2))]; // format("After remove: length(ind) = ``\n",length(ind)); MAP = MAP[ind,:]; }; function drawtitles(year,month,day,sec,N,F) global(longcoeff,lat1_map,lat2_map,long1_map,long2_map,monthstr) { hour = floor(sec/3600); minutes = floor(sec/60) - 60*hour; seconds = round(sec) - 3600*hour - 60*minutes; x = 0.5*(long1_map + long2_map); y = 0.5*(lat1_map + lat2_map); plot(longcoeff*#(x,x),#(y,y),"linetype",0, "toplabel",#("STARE ",where(N,"N",""),where(F,"F",""),where(!F && !N," (No data)","")), "subtitle",sformat("``-`` ``, UT ``:``:``",monthstr[month,:],sprintf("%02d",day), sprintf("%04d",where(year>1900,year,where(year>96,1900+year,2000+year))),sprintf("%02d",where(isint(hour),hour,0)), sprintf("%02d",minutes),sprintf("%02d",seconds)), "xlabel",sformat("Geogr. long * ``",longcoeff),"ylabel","Geogr. lat"); }; function drawmap() global(MAP,longcoeff, lat_malvik,long_malvik, lat_hankasalmi,long_hankasalmi, UseColorPrinting) { if (!isarray(MAP)) { format("Loading map file MAAMERI.HDF ... "); loadmapfile(); format("done.\n"); }; [n,two] = size(MAP); firsttime = 1; hold(on); for (i=1; i<=n; i++) { for (j=i; j<=n; j++) if (MAP[j,1] == -999.0) break; if (j > i) { lat = MAP[i:j-1,1]; long = -MAP[i:j-1,2]; plot(longcoeff*long,lat,"linecolor",where(UseColorPrinting,0,yellow)); firsttime = 0; }; i = j; }; annotate("point","x1",longcoeff*long_hankasalmi,"y1",lat_hankasalmi,"markercolor",red,"markertype",13,"markersize",3); annotate("point","x1",longcoeff*long_malvik,"y1",lat_malvik,"markercolor",red,"markertype",13,"markersize",3); hold(off); }; function u = Data(lat,long) global(lat1,lat2,long1,long2) { lat0 = 0.5*(lat1 + lat2); long0 = 0.5*(long1 + long2); u = exp(-((lat-lat0)/2)^2-((long-long0)^2)^2); }; function [x,y,z] = fromspherical(r,theta,phi) { x = r*sin(theta)*cos(phi); y = r*sin(theta)*sin(phi); z = r*cos(theta); }; function [r,theta,phi] = tospherical(x,y,z) { r = sqrt(x^2 + y^2 + z^2); phi = atan(y/x); theta = acos(z/r); }; function [lat,long,data] = initdata() global(lat1,lat2,long1,long2,longcoeff) { latstep = 0.2; longstep = 0.4; [lat,long] = grid(lat1:latstep:lat2, long1:longstep:long2); data = 0*lat; }; function [lat,long,data] = initdata_v() global(lat1,lat2,long1,long2,longcoeff) { latstep = 0.5; longstep = 1.0; [lat,long] = grid(lat1:latstep:lat2, long1:longstep:long2); data = 0*lat; }; function [data] = normalize(datacnt) { data*= 1/max(1,datacnt); }; function [lat,long;evec] = genbeam(lat0,long0,angle;range1) { R_E_km = 6378; if (isarray(range1)) range = range1 else { maxrange = 1100.0; rangestep = 10.0; range = rangestep:rangestep:maxrange; }; theta0 = (pi/180)*(90-lat0); phi0 = (pi/180)*long0; [x0,y0,z0] = fromspherical(R_E_km,theta0,phi0); etheta = #(cos(theta0)*cos(phi0),cos(theta0)*sin(phi0),-sin(theta0)); ephi = #(-sin(phi0),cos(phi0),0); evec = cos(angle)*(-etheta) + sin(angle)*ephi; x = x0 + range*evec[1]; y = y0 + range*evec[2]; z = z0 + range*evec[3]; [r,theta,phi] = tospherical(x,y,z); lat = 90 - (180/pi)*theta; long = (180/pi)*phi; }; function u = kernel(lat0,long0, nlat,nlong, lat,long) global(lat1,lat2) { dlat = lat - lat0; dlong = long - long0; wpar = 0.15; wperp = 0.4 + 0.6*(lat-lat1)/(lat2-lat1); // u = exp(-((nlat*dlat + nlong*dlong)/wpar)^2 - ((nlong*dlat - nlat*dlong)/wperp)^2); logu = -((lat-lat0)/wperp)^2 - ((long-long0)/wperp)^2; logu+= -((nlat*dlat + nlong*dlong)/wpar)^2; u = exp(logu); }; function [data,datacnt] = gatherdata1(u,range, lat,long, lat0,long0,angle) // Gather data for one beam. u and range must be vectors (of length 50, usually) // lat,long must be matrices of the same size as data,datacnt // lat0,long0,angle define the station and the beam direction { [latv,longv] = genbeam(lat0,long0,angle,range); n = length(latv); nlat = latv[1] - latv[2]; nlong = longv[1] - longv[2]; for (i=1; i<=n; i++) { if (i > 1) { nlat = latv[i-1] - latv[i]; nlong = cos(latv[i]*(pi/180))*(longv[i-1] - longv[i]); }; norm = 1/sqrt(nlat^2 + nlong^2); nlat*= norm; nlong*= norm; // format("nlat=``, nlong=``\n",nlat,nlong); w = kernel(latv[i],longv[i], nlat,nlong, lat,long); data+= w*u[i]; datacnt+= w; }; // format("gatherdata1: data=`` .. ``\n",min(data),max(data)); }; function [data_vx,data_vy,data_vz, datacnt_vx,datacnt_vy,datacnt_vz] = gatherdatav1(u,range, lat,long, lat0,long0,angle) { [latv,longv,evec] = genbeam(lat0,long0,angle,range); n = length(latv); nlat = latv[1] - latv[2]; nlong = longv[1] - longv[2]; for (i=1; i<=n; i++) { if (i > 1) { nlat = latv[i-1] - latv[i]; nlong = cos(latv[i]*(pi/180))*(longv[i-1] - longv[i]); }; norm = 1/sqrt(nlat^2 + nlong^2); nlat*= norm; nlong*= norm; w = kernel(latv[i],longv[i], nlat,nlong, lat,long); data_vx-= w*evec[1]*u[i]; data_vy-= w*evec[2]*u[i]; data_vz-= w*evec[3]*u[i]; datacnt_vx+= w; datacnt_vy+= w; datacnt_vz+= w; }; }; function [data,datacnt] = gatherdata(u,range, lat,long, lat0,long0,angles) { n = length(angles); for (i=1; i<=n; i++) [data,datacnt] = gatherdata1(u[i,:],range, lat,long, lat0,long0, angles[i]); }; function [data_vx,data_vy,data_vz, datacnt_vx,datacnt_vy,datacnt_vz] = gatherdatav(u,range, lat,long, lat0,long0,angles) { n = length(angles); for (i=1; i<=n; i++) [data_vx,data_vy,data_vz, datacnt_vx,datacnt_vy,datacnt_vz] = gatherdatav1(u[i,:],range,lat,long,lat0,long0,angles[i]); }; function drawdata(data) global(lat1,lat2,long1,long2,longcoeff,UseInterpolation,UseColorPrinting, PLOTTED_DATA, PLOTTED_LONG1,PLOTTED_LONG2,PLOTTED_LAT1,PLOTTED_LAT2) { PLOTTED_DATA = transpose(data); PLOTTED_LONG1 = long1; PLOTTED_LONG2 = long2; PLOTTED_LAT1 = lat1; PLOTTED_LAT2 = lat2; pcolor(transpose(data),"cmin",where(UseColorPrinting,-10,0),"cmax",30, "xmin",longcoeff*long1,"xmax",longcoeff*long2,"ymin",lat1,"ymax",lat2, "interpolate",where(UseInterpolation,1,3)); }; function drawbeam(lat0,long0,angle) global(longcoeff,UseColorPrinting) { [lat,long] = genbeam(lat0,long0,angle /*,495:15:1230*/); plot(longcoeff*long,lat,"linecolor",where(UseColorPrinting,0,yellow)); }; function drawbeams() global(lat_malvik,long_malvik,lat_hankasalmi,long_hankasalmi, angles_malvik,angles_hankasalmi) { hold(on); for (i=1; i<=length(angles_malvik); i++) drawbeam(lat_malvik,long_malvik,angles_malvik[i]); for (i=1; i<=length(angles_hankasalmi); i++) drawbeam(lat_hankasalmi,long_hankasalmi,angles_hankasalmi[i]); hold(off); }; function [range,sttime,entime,snr,vel, success] = readfile(year,month,day,NF) global(DataDir) { // Try to read in day-file year4 = year; if (year4 < 1900) { if (year4 > 96) year4+= 1900 else year4+= 2000; }; dir = DataDir; yearstr = sprintf("%04d",year4); monthstr = sprintf("%02d",month); daystr = sprintf("%02d",day); inputfile = #(yearstr,"/",monthstr,"/",yearstr,monthstr,daystr,".",where(NF=='n',"nor","fin")); // inputfile = #(sprintf("%04d",year4),sprintf("%02d",month),sprintf("%02d",day),".",where(NF=='n',"nor","fin")); inputfile_analyzed = #(inputfile,".analyzed.hdf"); inputfile_analyzed_gz = #(inputfile_analyzed,".gz"); delflag = 0; if (isfile(#(dir,"/",inputfile_analyzed_gz))) { if (!isfile(#(dir,"/",inputfile_analyzed))) { tmpdir = getenv("TMPDIR"); if (!isstring(tmpdir)) tmpdir = "/tmp"; inputfile_analyzed = #(tmpdir,"/",inputfile_analyzed[9:length(inputfile_analyzed)]); format("Uncompressing ``\n",inputfile_analyzed_gz); system(sformat("gunzip <`` >``",#(dir,"/",inputfile_analyzed_gz),inputfile_analyzed)); delflag = 1; } }; if (isfile(inputfile_analyzed)) { format("Reading analyzed file `` ...",inputfile_analyzed); range = import1(inputfile_analyzed,"RealArray:range"); if (min(range) < 0) { format("\n***** WARNING: min(range)=`` < 0\n",min(range)); format("***** Resetting range to 495:15:1230\n"); range = 495:15.0:1230; }; sttime = import1(inputfile_analyzed,"RealArray:sttime"); entime = import1(inputfile_analyzed,"RealArray:entime"); if (max(max(entime),max(sttime)) > 1.01*24*3600) { format("\n***** WARNING: sttime,entime in analyzed file seems to be counted\n"); format(" from the beginning of YEAR, not DAY as it should.\n"); format(" Doing st/entime=st/entime-min(st/entime); this will give incorrect\n"); format(" result if the day starts with a data gap.\n"); sttime = sttime - min(sttime); entime = entime - min(entime); }; snr = import1(inputfile_analyzed,"RealArray:snr"); vel = import1(inputfile_analyzed,"RealArray:vel"); snr = max(1,snr); success = isarray(entime) && isarray(range) && isarray(sttime) && isarray(snr) && isarray(vel); if (!success) format("no success\n") else format("ok\n"); if (delflag) remove(inputfile_analyzed); // IF the file was produced by uncompression, remove it now // This is A BIT dangerous, but it SEEMS to work, and it is ANALYZED data anyway... } else { success = 0; }; }; // for plotit() and drawvectors(): set doplot=0 if you just want the PLOTTED_VLONG etc. out of it function drawvectors(vx,vy,vz, lat,long) global(longcoeff,lat1,lat2,long1,long2, PLOTTED_VLONG, PLOTTED_VLAT, PLOTTED_VGRID_LONG, PLOTTED_VGRID_LAT, doplot, PlotEField) { [N,M] = size(lat); vlong = 0*lat; vlat = vlong; phi = (pi/180)*long; theta = (pi/180)*(90-lat); ephi_x = -sin(phi); ephi_y = cos(phi); ephi_z = 0; etheta_x = cos(theta)*cos(phi); etheta_y = cos(theta)*sin(phi); etheta_z = -sin(theta); vlong = ephi_x*vx + ephi_y*vy + ephi_z*vz; vlat = -(etheta_x*vx + etheta_y*vy + etheta_z*vz); if (PlotEField) { format("Plotting electric field.\n"); format("One degree corresponds to 50 mV/m.\n"); Babs = 50000e-9; er_x = sin(theta)*cos(phi); er_y = sin(theta)*sin(phi); er_z = cos(theta); // B = - |B| er Bx = -Babs*er_x; By = -Babs*er_y; Bz = -Babs*er_z; // E = B x v Ex = By*vz - Bz*vy; Ey = Bz*vx - Bx*vz; Ez = Bx*vy - By*vx; Elong = ephi_x*Ex + ephi_y*Ey + ephi_z*Ez; Elat = -(etheta_x*Ex + etheta_y*Ey + etheta_z*Ez); PLOTTED_VLONG = Elong; PLOTTED_VLAT = Elat; vscale = 20; // one degree corresponds to 0.05 V/m electric field } else { format("Plotting electron velocity field.\n"); format("One degree corresponds to 1 km/s.\n"); PLOTTED_VLONG = vlong; PLOTTED_VLAT = vlat; vscale = 0.001; // one degree corresponds to 1000 m/s velocity }; PLOTTED_VGRID_LONG = long; PLOTTED_VGRID_LAT = lat; if (doplot) { hold(on); plot(longcoeff*#(floor(long1-1),ceil(long2+1)),#(floor(lat1-1),ceil(lat2+1)),"linestyle",0); vplot(longcoeff*long,lat,/*longcoeff* */PLOTTED_VLONG,PLOTTED_VLAT,"vscale",vscale); hold(off); } }; function beampermF = FinnishBeamSwap(year,month,day) { if (year < 1998 || (year==1998 && (month < 6 || month==6 && day<17))) { format("Swapping Finnish beams 3 and 6\n"); beampermF = #(1,2,6,4,5,3,7,8); } else beampermF = 1:8; }; function plotit(rangeN,rangeF, sttimeN,entimeN, sttimeF,entimeF, snrN,snrF,velN,velF) global(lat_malvik,long_malvik,lat_hankasalmi,long_hankasalmi,angles_malvik,angles_hankasalmi, successF,successN, data_snrN,data_snrF,vx,vy,vz, year,month,day,time_wanted,integration_time, NorwegianSign,FinnishSign, DisplayMap,DisplayBeams,DisplayFinnish,DisplayNorwegian,DisplaySNR,DisplayVectors,DisplayBothVectors,doplot, DBLIMIT,UseStrictDBLimit) { [nbeams,nranges,ntimes] = size(snrN); beampermF = FinnishBeamSwap(year,month,day); [lat,long,data_snrN] = initdata(); data_snrF = data_snrN; datacnt_snrN = data_snrN; datacnt_snrF = data_snrF; // time_wanted contains time from start of UT-day in seconds // now find the corresponding timeindF,timeindN for both radars [mindtF,timeindF] = min(abs(time_wanted - 0.5*(sttimeF+entimeF))); [mindtN,timeindN] = min(abs(time_wanted - 0.5*(sttimeN+entimeN))); gapF = 0; gapN = 0; if (mindtF > 0.5*integration_time) { format("- Data gap in Finnish data, mindtF=``\n",mindtF); gapF = 1; }; if (mindtN > 0.5*integration_time) { format("- Data gap in Norwegian data, mindtN=``\n",mindtN); gapN = 1; }; // format("- timeindF = ``, timeindN = ``\n",timeindF,timeindN); if (DisplaySNR) { if (DisplayNorwegian && successN && !gapN) { // format("Gathering snrN..\n"); [data_snrN,datacnt_snrN] = gatherdata(snrN[:,:,timeindN],rangeN, lat,long, lat_malvik,long_malvik,angles_malvik); }; if (DisplayFinnish && successF && !gapF) { // format("Gathering snrF..\n"); [data_snrF,datacnt_snrF] = gatherdata(snrF[beampermF,:,timeindF],rangeF, lat,long, lat_hankasalmi,long_hankasalmi,angles_hankasalmi); // format("snrF: `` .. ``\n",min(snrF[:,:,timeindF]),max(snrF[:,:,timeindF])); // format("max(data_snrF) = ``\n",max(data_snrF)); // format("max(datacnt_snrF) = ``\n",max(datacnt_snrF)); }; [data_snrN] = normalize(datacnt_snrN); [data_snrF] = normalize(datacnt_snrF); // format("data_snrN: `` .. ``\n",min(data_snrN),max(data_snrN)); // format("data_snrF: `` .. ``\n",min(data_snrF),max(data_snrF)); }; if (DisplayVectors) { [lat_v,long_v,data_vxN] = initdata_v(); data_vyN = data_vxN; data_vzN = data_vyN; datacnt_vxN = data_vxN; datacnt_vyN = data_vyN; datacnt_vzN = data_vzN; data_snrN_2 = data_vxN; data_snrF_2 = data_vxN; datacnt_snrN_2 = data_vxN; datacnt_snrF_2 = data_vxN; data_vxF = data_vxN; data_vyF = data_vxN; data_vzF = data_vxN; datacnt_vxF = data_vxN; datacnt_vyF = data_vxN; datacnt_vzF = data_vxN; hasNdata = 0; hasFdata = 0; if ((DisplayBothVectors || DisplayNorwegian) && successN && !gapN) { [data_vxN,data_vyN,data_vzN, datacnt_vxN,datacnt_vyN,datacnt_vzN] = gatherdatav(NorwegianSign*velN[:,:,timeindN],rangeN,lat_v,long_v, lat_malvik,long_malvik,angles_malvik); [data_vxN] = normalize(datacnt_vxN); [data_vyN] = normalize(datacnt_vyN); [data_vzN] = normalize(datacnt_vzN); hasNdata = 1; }; if ((DisplayBothVectors || DisplayFinnish) && successF && !gapF) { [data_vxF,data_vyF,data_vzF, datacnt_vxF,datacnt_vyF,datacnt_vzF] = gatherdatav(FinnishSign*velF[beampermF,:,timeindF],rangeF,lat_v,long_v, lat_hankasalmi,long_hankasalmi,angles_hankasalmi); [data_vxF] = normalize(datacnt_vxF); [data_vyF] = normalize(datacnt_vyF); [data_vzF] = normalize(datacnt_vzF); hasFdata = 1; }; if (hasNdata) { if (hasFdata) { // Both N and F data (normal case) N2 = data_vxN^2 + data_vyN^2 + data_vzN^2; F2 = data_vxF^2 + data_vyF^2 + data_vzF^2; dot = data_vxN*data_vxF + data_vyN*data_vyF + data_vzN*data_vzF; badind = find(N2*F2 == dot^2); // these will produce division by zero, remove them soon... coeffN = (N2 - dot)*F2/(N2*F2 - dot^2); coeffF = (F2 - dot)*N2/(N2*F2 - dot^2); data_vx = coeffN*data_vxN + coeffF*data_vxF; data_vy = coeffN*data_vyN + coeffF*data_vyF; data_vz = coeffN*data_vzN + coeffF*data_vzF; data_vx[badind] = 0; data_vy[badind] = 0; data_vz[badind] = 0; } else { // N data but no F data data_vx = data_vxN; data_vy = data_vyN; data_vz = data_vzN; } } else { if (hasFdata) { // F data but no N data data_vx = data_vxF; data_vy = data_vyF; data_vz = data_vzF; } else { // Neither N nor F data data_vx = 0*data_vxN; data_vy = 0*data_vyN; data_vz = 0*data_vzN; }; }; // For too small SNR values, set the corresponding vectors to zero // For this, we need to accumulate SNR again since vector grid and SNR grid have different sizes... if (successN && !gapN) [data_snrN_2, datacnt_snrN_2] = gatherdata(snrN[:,:,timeindN],rangeN, lat_v,long_v, lat_malvik,long_malvik,angles_malvik); if (successF && !gapF) [data_snrF_2, datacnt_snrF_2] = gatherdata(snrF[:,:,timeindF],rangeF, lat_v,long_v, lat_hankasalmi,long_hankasalmi,angles_hankasalmi); [data_snrN_2] = normalize(datacnt_snrN_2); [data_snrF_2] = normalize(datacnt_snrF_2); // Now zero [n,m] = size(data_vx); zerocnt = 0; if (successN && !gapN && successF && !gapF) { if (UseStrictDBLimit) data_crit = min(data_snrN_2, data_snrF_2) else data_crit = 0.5*(data_snrN_2 + data_snrF_2); } else data_crit = data_snrN_2 + data_snrF_2; for (i=1; i<=n; i++) for (j=1; j<=m; j++) if (data_crit[i,j] < DBLIMIT) { data_vx[i,j] = 0; data_vy[i,j] = 0; data_vz[i,j] = 0; zerocnt++; }; // format("Zeroed `` vectors\n",zerocnt); }; if (doplot) { hold(on); drawtitles(year,month,day,time_wanted /*sttimeF[timeind]*/ /*- sttimeF[1]*/, DisplayNorwegian && successN && !gapN, DisplayFinnish && successF && !gapF); if (DisplayMap) drawmap(); if (successN && !gapN && successF && !gapF) // then we have both radars, plot the average snr_plotted = 0.5*(data_snrN + data_snrF) else snr_plotted = data_snrN + data_snrF; // otherwise we take the sum rather than average if (DisplaySNR && max(snr_plotted) > 0) drawdata(snr_plotted); if (DisplayBeams) drawbeams(); if (DisplayVectors) drawvectors(data_vx,data_vy,data_vz,lat_v,long_v); hold(off); } else { drawvectors(data_vx,data_vy,data_vz,lat_v,long_v); }; }; function success=ReadDataFiles() global { [rangeN,sttimeN,entimeN,snrN,velN, successN] = readfile(year,month,day,'n'); [rangeF,sttimeF,entimeF,snrF,velF, successF] = readfile(year,month,day,'f'); // if (successN || successF) integration_time = sttimeF[2] - sttimeF[1]; if (successN) { integration_time = sttimeN[2] - sttimeN[1]; } else if (successF) { integration_time = sttimeF[2] - sttimeF[1]; }; // maxtimeindN = 0; maxtimeindF = 0; // if (successN) maxtimeindN = length(sttimeN); // if (successF) maxtimeindF = length(sttimeF); // maxtimeind = max(maxtimeindF,maxtimeindN); if (!successN && !successF) { format("*** Neither Finnish nor Norwegian data could be read\n"); } else { if (!successF) { format("Warning: No Finnish data\n"); rangeF = rangeN; snrF = 0*snrN; velF = 0*velN; sttimeF = sttimeN; entimeF = entimeN; }; if (!successN) { format("Warning: No Norwegian data\n"); rangeN = rangeF; snrN = 0*snrF; velN = 0*velF; sttimeN = sttimeF; entimeN = entimeF; }; }; success = successN || successF; }; function UpdatePlot() global { if (success) plotit(rangeN,rangeF, sttimeN,entimeN, sttimeF,entimeF, snrN,snrF,velN,velF) else format("*** No date loaded, cannot plot\n"); }; function UpdateColorMode() global { if (UseColorPrinting) plotopt("-colorps -bg white -fg black") else plotopt(""); }; function RTIplot(snr,sttime,entime,toplabel,year,month,day) global(monthstr) { beam = 5; [c8,nranges,ntimes1] = size(snr); tAhour = 0; tBhour = 24; // tAhour = 21; // tBhour = 24; if (tBhour-tAhour < 12) interleave = 1 else interleave = 3*3; timeres = 20; data = zeros(round((tBhour-tAhour)*3600/(timeres*interleave)),nranges); [ntimes,dummy] = size(data); for (i=1; i<=ntimes; i++) { tsec = tAhour*3600 + (i-1)*timeres*interleave; [mindt,minpos] = min(abs(tsec-0.5*(sttime+entime))); if (mindt <= 0.5*timeres) data[i,:] = snr[beam,:,minpos]; }; pcolor(data,"toplabel",toplabel,"equalscale","false","xmin",tAhour,"xmax",tBhour,"xlabel","UT (hour)","ylabel","", "cmin",0,"cmax",ceil(max(data)), "comment",sformat("RTI plot for beam ``",beam), "subtitle",sformat("``-`` ``",monthstr[month,:],sprintf("%02d",day), sprintf("%04d",where(year > 1900, year, where(year>96,1900+year,2000+year))))); }; function RTIplots() global(snrF,snrN,sttimeF,sttimeN,entimeF,entimeN,year,month,day) { oldmode = holdmode(); holdmode(stacking); hold(on); RTIplot(snrF,sttimeF,entimeF,"STARE FINLAND",year,month,day); RTIplot(snrN,sttimeN,entimeN,"STARE NORWAY",year,month,day); hold(off); holdmode(oldmode); }; function dump_LTV(duration,integration,ltv_dir) local(dur,finaltime,integration_time,cnt,Ntimes,orig_time_wanted) { // Making sure that the filename is proper year_str=sformat("``",year); if (month < 10 ) { month_str=sformat("0``",month); } else { month_str=sformat("``",month); }; if (day < 10) { day_str=sformat("0``",day); } else { day_str=sformat("0``",day); }; datum=#(year_str,month_str,day_str,"LTV.hdf"); // a full day's worth of data dur = duration; integration_time = integration; finaltime = time_wanted + 60*dur; doplot = 0; // Temporarily disable plotting, just gather PLOTTED_VLONG,PLOTTED_VLAT,PLOTTED_VGRID_LONG,PLOTTED_VGRID_LAT cnt = 0; // count number of times first orig_time_wanted = time_wanted; repeat time_wanted = limit(time_wanted+integration_time,0,24*3600); cnt++; until time_wanted > finaltime || time_wanted >= 24*3600; Ntimes = cnt; time_wanted = orig_time_wanted; // format("`` time steps\n",Ntimes); cnt = 0; repeat UpdatePlot(); if (cnt == 0) { [c11,c14] = size(PLOTTED_VLAT); LTV_Vnorth = zeros(Ntimes,c11,c14); LTV_Veast = LTV_Vnorth; LTV_time = zeros(Ntimes); }; // format("t = `` ready\n",time_wanted); cnt++; LTV_Vnorth[cnt,:,:] = PLOTTED_VLAT; LTV_Veast[cnt,:,:] = PLOTTED_VLONG; LTV_time[cnt] = time_wanted; time_wanted = limit(time_wanted+integration_time,0,24*3600); until time_wanted > finaltime || time_wanted >= 24*3600; doplot = 1; // ltv_dir is the directory where the gzipped data should be saved cd(ltv_dir); save(datum,"PLOTTED_VGRID_LAT","PLOTTED_VGRID_LONG","LTV_Vnorth","LTV_Veast","LTV_time"); compress_ltv=sformat("gzip ``",datum); system(compress_ltv); // save("LTV.hdf","PLOTTED_VGRID_LAT","PLOTTED_VGRID_LONG","LTV_Vnorth","LTV_Veast","LTV_time"); }; function RTIBatchAnalysisOperation(analyzedfile) global(year,month,day,DataDir) { /* * tela -be 'analyzedfile="1997/12/19971202.fin.analyzed.hdf.gz"' stare.t * This command writes /proj/STARE/RTI/1997/12/19971202.RTI.gif. */ year = evalexpr(analyzedfile[1:4]); month = evalexpr(analyzedfile[6:7]); day = evalexpr(analyzedfile[15:16]); format("year=``, month=``, day=``\n",year,month,day); cd(DataDir); success = ReadDataFiles(); if (!success) return; plotopt("-bg white -fg black -colorps -o dataplot.ps -noxplot -print -printcmd echo -Pnull"); remove("dataplot.ps"); RTIplots(); // Wait until dataplot.ps appears repeat if (isfile("dataplot.ps")) break; pause(1); until 0; // Then write 5 seconds so that the file gets written properly ... pause(5); // Change the string "Times-" to "Charter-" in dataplot.ps. // The Charter font is the only TrueType font known to Ghostscript. system("sed -e 's/Times-/Charter-/g' x.ps; mv x.ps dataplot.ps"); system("echo quit | gs -sDEVICE=ppm -sOutputFile=dataplot.ppm -r68 dataplot.ps"); outputGIFfile = #("../RTI/",analyzedfile[1:7],"/",analyzedfile[9:16],".RTI.gif"); system(sformat("ppmtogif dataplot.ppm >``; chmod og+r ``; rm dataplot.ppm", outputGIFfile,outputGIFfile,outputGIFfile)); }; function MakeAnimation() local(dur,finaltime,display,current_display,integration_time) { format("Making animation from current time forward.\n"); dur = input("Give duration of movie in minutes"); integration_time = input("Give time interval between frames in seconds (e.g. 20)"); format("The animation will contain ~ `` frames.\n",round(dur*60/integration_time)); system("rm -f stare-anim-????.gif"); finaltime = time_wanted + 60*dur; current_display = getenv("DISPLAY"); display = current_display; if (!isstring(current_display)) current_display = ""; ans = menu("Select X display", where(streq(current_display,""),"(empty string)",current_display), "esim.fmi.fi:2", "nabla.fmi.fi:1", "Other"); if (ans == 1) { display = current_display; } else if (ans == 2) { display = "esim.fmi.fi:2"; } else if (ans == 3) { display = "nabla.fmi.fi:1"; } else if (ans == 4) { display = input_string("Give X display name:"); }; set_alt_display(display); // plotopt("-bg white -fg black -colorps -gifsave -o stare-anim-#.gif -printcmd echo -Pnull"); plotopt(sformat("-display `` -bg white -fg black -colorps -gifsave -o stare-anim-#.gif",display)); closefig(1); cnt = 0; repeat UpdatePlot(); wait_time = 0; repeat pause(1); wait_time++; until isfile(sprintf("stare-anim-%.4d.gif",cnt)) || (wait_time > 20); pause(1); wait_time++; cnt++; // pause(4); format("t = `` ready, waited for `` seconds\n",time_wanted,wait_time); time_wanted = limit(time_wanted+integration_time,0,24*3600); until time_wanted > finaltime || time_wanted >= 24*3600; s = "makemovie -c qt_anim -r 1 -o stareanim.qt -f qt stare-anim-????.gif"; format(s); system(s); system("rm stare-anim-????.gif"); format("Output file is stareanim.qt.\n"); reset_alt_display(); UpdateColorMode(); closefig(1); UpdatePlot(); }; function LTVplot() local(dur,finaltime,integration_time,cnt,Ntimes,orig_time_wanted) { format("Making LTV plot from current time forward.\n"); dur = input("Give duration of event in minutes"); integration_time = input("Give wanted time interval between frames in seconds (e.g. 20)"); finaltime = time_wanted + 60*dur; doplot = 0; // Temporarily disable plotting, just gather PLOTTED_VLONG,PLOTTED_VLAT,PLOTTED_VGRID_LONG,PLOTTED_VGRID_LAT cnt = 0; // count number of times first orig_time_wanted = time_wanted; repeat time_wanted = limit(time_wanted+integration_time,0,24*3600); cnt++; until time_wanted > finaltime || time_wanted >= 24*3600; Ntimes = cnt; time_wanted = orig_time_wanted; format("`` time steps\n",Ntimes); cnt = 0; repeat UpdatePlot(); if (cnt == 0) { [c11,c14] = size(PLOTTED_VLAT); LTV_Vnorth = zeros(Ntimes,c11,c14); LTV_Veast = LTV_Vnorth; LTV_time = zeros(Ntimes); }; format("t = `` ready\n",time_wanted); cnt++; LTV_Vnorth[cnt,:,:] = PLOTTED_VLAT; LTV_Veast[cnt,:,:] = PLOTTED_VLONG; LTV_time[cnt] = time_wanted; time_wanted = limit(time_wanted+integration_time,0,24*3600); until time_wanted > finaltime || time_wanted >= 24*3600; doplot = 1; save("LTV.hdf","PLOTTED_VGRID_LAT","PLOTTED_VGRID_LONG","LTV_Vnorth","LTV_Veast","LTV_time"); }; function GrabLastToHdfFile() local(fn,hh,minutes,ss) { hh = floor(time_wanted/3600); minutes = floor(time_wanted/60 - 60*hh); ss = round(time_wanted - 60*minutes - 3600*hh); fn = #(sprintf("%04d",year), sprintf("%02d",month), sprintf("%02d",day), ".", sprintf("%02d",hh), sprintf("%02d",minutes), sprintf("%02d",ss), ".vars.hdf"); ans = input(sformat("Give output filename [``] ",fn)); if (length(ans) > 0) fn = ans; save(fn, "PLOTTED_DATA", "PLOTTED_LAT1","PLOTTED_LAT2", "PLOTTED_LONG1","PLOTTED_LONG2", "PLOTTED_VGRID_LAT","PLOTTED_VGRID_LONG", "PLOTTED_VLAT","PLOTTED_VLONG"); format("Wrote \"``\"\n",fn); }; function GrabLastToAsciiFile() local(fn,hh,minutes,ss,fp,N,i,ans) { hh = floor(time_wanted/3600); minutes = floor(time_wanted/60 - 60*hh); ss = round(time_wanted - 60*minutes - 3600*hh); fn = #(sprintf("%04d",year), sprintf("%02d",month), sprintf("%02d",day), ".", sprintf("%02d",hh), sprintf("%02d",minutes), sprintf("%02d",ss), ".dat"); ans = input(sformat("Give output filename [``] ",fn)); if (length(ans) > 0) fn = ans; fp = fopen(fn,"w"); N = length(PLOTTED_VGRID_LAT); if (PlotEField) scaling = 1000 else scaling = 1; for (i=1; i<=N; i++) if (PLOTTED_VLAT[i] != 0 || PLOTTED_VLONG[i] != 0) fformat(fp,"`` `` `` ``\n",PLOTTED_VGRID_LAT[i],PLOTTED_VGRID_LONG[i],scaling*PLOTTED_VLAT[i],scaling*PLOTTED_VLONG[i]); fclose(fp); format("Wrote `` field to ASCII file \"``\"\n",where(PlotEField,"electric","velocity"),fn); if (PlotEField) format("Format of file: lat long Enorth(mV/m) Eeast(mV/m)\n") else format("Format of file: lat long Vnorth(m/s) Veast(m/s)\n") }; function PrintDataAt(atlat,atlon) local(latstep,longstep,data_snrF,data_snrN,datacnt_snrF,datacnt_snrN, gapF,gapN) { lat = zeros(1,1); long = zeros(1,1); lat[1,1] = atlat; long[1,1] = atlon; data_snrN = 0*lat; beampermF = FinnishBeamSwap(year,month,day); data_snrF = data_snrN; datacnt_snrN = data_snrN; datacnt_snrF = data_snrF; data_vxN = data_snrN; data_vyN = data_vxN; data_vzN = data_vyN; datacnt_vxN = data_vxN; datacnt_vyN = data_vyN; datacnt_vzN = data_vzN; data_snrN_2 = data_vxN; data_snrF_2 = data_vxN; datacnt_snrN_2 = data_vxN; datacnt_snrF_2 = data_vxN; data_vxF = data_vxN; data_vyF = data_vxN; data_vzF = data_vxN; datacnt_vxF = data_vxN; datacnt_vyF = data_vxN; datacnt_vzF = data_vxN; // time_wanted contains time from start of UT-day in seconds // now find the corresponding timeindF,timeindN for both radars [mindtF,timeindF] = min(abs(time_wanted - 0.5*(sttimeF+entimeF))); [mindtN,timeindN] = min(abs(time_wanted - 0.5*(sttimeN+entimeN))); gapF = 0; gapN = 0; if (mindtF > 0.5*integration_time) { format("- Data gap in Finnish data, mindtF=``\n",mindtF); gapF = 1; }; if (mindtN > 0.5*integration_time) { format("- Data gap in Norwegian data, mindtN=``\n",mindtN); gapN = 1; }; hasNdata = 0; hasFdata = 0; if (!gapN) { [data_snrN,datacnt_snrN] = gatherdata(snrN[:,:,timeindN],rangeN, lat,long, lat_malvik,long_malvik,angles_malvik); [data_vxN,data_vyN,data_vzN, datacnt_vxN,datacnt_vyN,datacnt_vzN] = gatherdatav(NorwegianSign*velN[:,:,timeindN],rangeN,lat,long, lat_malvik,long_malvik,angles_malvik); [data_vxN] = normalize(datacnt_vxN); [data_vyN] = normalize(datacnt_vyN); [data_vzN] = normalize(datacnt_vzN); hasNdata = 1; }; if (!gapF) { [data_snrF,datacnt_snrF] = gatherdata(snrF[beampermF,:,timeindF],rangeF, lat,long, lat_hankasalmi,long_hankasalmi,angles_hankasalmi); [data_vxF,data_vyF,data_vzF, datacnt_vxF,datacnt_vyF,datacnt_vzF] = gatherdatav(FinnishSign*velF[beampermF,:,timeindF],rangeF,lat,long, lat_hankasalmi,long_hankasalmi,angles_hankasalmi); [data_vxF] = normalize(datacnt_vxF); [data_vyF] = normalize(datacnt_vyF); [data_vzF] = normalize(datacnt_vzF); hasFdata = 1; }; [data_snrN] = normalize(datacnt_snrN); [data_snrF] = normalize(datacnt_snrF); if (hasNdata) { if (hasFdata) { // Both N and F data (normal case) N2 = data_vxN^2 + data_vyN^2 + data_vzN^2; F2 = data_vxF^2 + data_vyF^2 + data_vzF^2; dot = data_vxN*data_vxF + data_vyN*data_vyF + data_vzN*data_vzF; badind = find(N2*F2 == dot^2); // these will produce division by zero, remove them soon... coeffN = (N2 - dot)*F2/(N2*F2 - dot^2); coeffF = (F2 - dot)*N2/(N2*F2 - dot^2); data_vx = coeffN*data_vxN + coeffF*data_vxF; data_vy = coeffN*data_vyN + coeffF*data_vyF; data_vz = coeffN*data_vzN + coeffF*data_vzF; data_vx[badind] = 0; data_vy[badind] = 0; data_vz[badind] = 0; } else { // N data but no F data data_vx = data_vxN; data_vy = data_vyN; data_vz = data_vzN; } } else { if (hasFdata) { // F data but no N data data_vx = data_vxF; data_vy = data_vyF; data_vz = data_vzF; } else { // Neither N nor F data data_vx = 0*data_vxN; data_vy = 0*data_vyN; data_vz = 0*data_vzN; }; }; vx = data_vx[1,1]; vy = data_vy[1,1]; vz = data_vz[1,1]; phi = (pi/180)*atlon; theta = (pi/180)*(90-atlat); ephi_x = -sin(phi); ephi_y = cos(phi); ephi_z = 0; etheta_x = cos(theta)*cos(phi); etheta_y = cos(theta)*sin(phi); etheta_z = -sin(theta); vlong = ephi_x*vx + ephi_y*vy + ephi_z*vz; vlat = -(etheta_x*vx + etheta_y*vy + etheta_z*vz); Babs = 50000e-9; er_x = sin(theta)*cos(phi); er_y = sin(theta)*sin(phi); er_z = cos(theta); // B = - |B| er Bx = -Babs*er_x; By = -Babs*er_y; Bz = -Babs*er_z; // E = B x v Ex = By*vz - Bz*vy; Ey = Bz*vx - Bx*vz; Ez = Bx*vy - By*vx; Elong = ephi_x*Ex + ephi_y*Ey + ephi_z*Ez; Elat = -(etheta_x*Ex + etheta_y*Ey + etheta_z*Ez); format("SNR_N = `` dB, SNR_F = `` dB\n",data_snrN[1,1],data_snrF[1,1]); format("vx_N = ``, vy_N = ``, vz_N = `` [m/s]\n",data_vxN[1,1],data_vyN[1,1],data_vzN[1,1]); format("vx_F = ``, vy_F = ``, vz_F = `` [m/s]\n",data_vxF[1,1],data_vyF[1,1],data_vzF[1,1]); format("vx = ``, vy = ``, vz = `` [m/s]\n",data_vx[1,1],data_vy[1,1],data_vz[1,1]); format("v_North = ``, v_West = `` [m/s]\n",vlat,-vlong); format("E_North = ``, E_West = `` [mV/m]\n",1000*Elat,-1000*Elong); };