; file: makedopsumcubes.pro for 2005-07-13-spot/Cubes ; last: May 20 2012 Rob Rutten @ Hollum ; note: uses assoc to combine the two types of substring-specified cubes ; ##### specify the input = output directory indir='../Cubes' ; @ adapt? ; no slash at end ; ##### specify the pair(s) of input cubes, blue wing first inwavident=['ha-0500','ha+0500'] ; ##### specify names for the output cubes; same nr characters, dopp, wing outwavident=['ha-dopp','ha-wing'] ; ========== should be parameterfree from here on ============= ; get all pairs of file names containing the inwavidents infile1=file_search(indir+'/*'+inwavident[0]+'*') infile2=file_search(indir+'/*'+inwavident[1]+'*') ; get number of pairs sizein1list=size(infile1) nfiles1=sizein1list(1) sizein2list=size(infile2) nfiles2=sizein2list(1) if (nfiles2 ne nfiles1) then print,' ===== error: input files not paired' ; read, process, write pairs of cubes for ifile=0, nfiles1-1 do begin ; get infilenames infilename1=infile1(ifile) infilename2=infile2(ifile) ; define output filenames pos=strpos(infilename1,inwavident[0]) outfilename1=infilename1 strput,outfilename1,outwavident[0],pos outfilename2=infilename2 strput,outfilename2,outwavident[1],pos ; get cube geometry and file datatype getdotcubegeom,infilename1,nx,ny,nt,datatype getdotcubegeom,infilename2,nx2,ny2,nt2,datatype2 if (nx2 ne nx or ny2 ne ny or nt2 ne nt or datatype2 ne datatype) then $ print,' === error: input cubes not the same geometry or data type' ; open input file 1 for assoc get_lun, inunit1 openr,inunit1,infilename1 if (datatype eq 'float') then inarr1=fltarr(nx,ny) if (datatype eq 'integer') then inarr1=intarr(nx,ny) if (datatype eq 'byte') then inarr1=bytarr(nx,ny) inassoc1 = assoc(inunit1,inarr1) ; open input file 2 for assoc get_lun, inunit2 openr,inunit2,infilename2 if (datatype eq 'float') then inarr2=fltarr(nx,ny) if (datatype eq 'integer') then inarr2=intarr(nx,ny) if (datatype eq 'byte') then inarr2=bytarr(nx,ny) inassoc2 = assoc(inunit2,inarr2) ; open output file 1 for assoc get_lun, outunit1 openw,outunit1,outfilename1 if (datatype eq 'float') then outarr1=fltarr(nx,ny) if (datatype eq 'integer') then outarr1=intarr(nx,ny) if (datatype eq 'byte') then outarr1=bytarr(nx,ny) outassoc1 = assoc(outunit1,outarr1) ; open output file 2 for assoc get_lun, outunit2 openw,outunit2,outfilename2 if (datatype eq 'float') then outarr2=fltarr(nx,ny) if (datatype eq 'integer') then outarr2=intarr(nx,ny) if (datatype eq 'byte') then outarr2=bytarr(nx,ny) outassoc2 = assoc(outunit2,outarr2) ; first pass to find extrema doppim=fltarr(nx,ny) maxdopp=-1000. mindopp=+1000. maxwing=0. for it=0,nt-1 do begin a=float(inassoc1(it)) b=float(inassoc2(it)) sumim=a+b maxsum=max(sumim) if (maxsum gt maxwing) then maxwing=maxsum doppim(*,*)=(a-b)/sumim minim=min(doppim) maxim=max(doppim) if (minim lt mindopp) then mindopp=minim if (maxim gt maxdopp) then maxdopp=maxim endfor ;; print,' mindopp, maxdopp =',mindopp,maxdopp if (-mindopp gt maxdopp) then maxdopp=-mindopp ; second pass to make and write Doppler and wingsum cubes for it=0,nt-1 do begin a=float(inassoc1(it)) b=float(inassoc2(it)) sumim=a+b doppim(*,*)=(a-b)/sumim outassoc1[it]=fix(doppim*15000./maxdopp) ; fill half integer range outassoc2[it]=fix(sumim+15000./maxwing) endfor free_lun,inunit1 free_lun,inunit2 free_lun,outunit1 free_lun,outunit2 ; show completion print,' wrote '+outfilename1 print,' wrote '+outfilename2 endfor ; end loop over pairs end