FUNCTION MEAN, A B = TOTAL(A)/N_ELEMENTS(A) RETURN, B & END ;The executing program is xyoff_tex.pro. It performs a 2D ;cross-correlation, yielding the lag of maximum correlation, the ;coefficient of maximum correlation and the coefficient at lag 0. In ;practice, I usually used lag 0. Within the limits of floating point ;computation, the cross correlation at 0 shift of an image and itself is ;1.0 and the shifts are 0.0 and 0.0 as would be expected. xyoff_tex.pro ;calls max_pos.pro to find the maximum position within an array...the ;cross-correlation function in this case. rca_2D.pro is a setup and ;calling program for xyoff_tex.pro. You may not need it. ; ;If in running these, you find that they call other programs that are ;not included, let Tim know and he will supply you with them. Likewise, ;if you have any other problems. ;xyoff_tex.pro: REMOVE THIS LINE ;+ ; NAME: ; XYOFF ; PURPOSE: ; Calculates the correlation function between two images using a box of ; size (bx,by) surrounding a point in the image. By default this point ; is the center of the image, (nx/2,ny/2). The routine calculates the ; offsets by finding the maximum of the correlation fucntion. ; CATEGORY: ; IMAGE PROCESSING ; CALLING SEQUENCE: ; s = xyoff(a,b,bx,by,hx,hy,xx,yy,amax,fna,wdow=wdow,mask=mask,/sobel,ccf=ccf ; one=one,xy,amax) ; INPUTS: ; a 2D reference image. ; b 2D image that you want to track with respect to a ; bx size of the box in the x-dimension that you wish ; to use to generate the cross correlation function ; by size of the box in the y-dimension that you wish ; to use to generate the cross correlation function. ; hx x position in the image that corresponds to the center ; of the box. (Optional) ; hy y position in the image that corresponds to the center ; of the box. (Optional) ; fna fourier transform of the 256x256 (Optional) ; subtended image from the center of a. ; KEYWORD PARAMETERS: ; wdow=wdow 2-d apodization window. (Optional; default=1) ; mask=mask 2-d mask to be applied in fourier space. (Optional;default=1) ; /SOBEL = If specified, then SOBEL keyword causes this routine ; to track on the gradient of the image instead of the ; image itself. (It uses a Sobel edge enhancement operator) ; OUTPUTS: ; s = A vector containing the x and y offsets that should be applied ; to an image to maximize the cross correlation. ; s(0) = The x offset. ; s(1) = The y offset. ; s(2) = max value of the correlation ; s(3) = cross correlation at 0 shift ; s(4) = maximum of the real part of absolute value of cross correlation ; NOTES: ; The executing program is xyoff.pro. It performs a 2D ; cross-correlation, yielding the lag of maximum correlation, the ; coefficient of maximum correlation and the coefficient at lag 0. In ; practice, I usually used lag 0. Within the limits of floating point ; computation, the cross correlation at 0 shift of an image and itself is ; 1.0 and the shifts are 0.0 and 0.0 as would be expected. xyoff_tex.pro ; calls max_pos.pro to find the maximum position within an array...the ; cross-correlation function in this case. ; MODIFICATION HISTORY: ; H. Cohl, 25 Mar, 1992 --- Made documentation more understandable. ; H. Cohl, 13 Feb, 1992 --- Made it so that you can enter in hx or hy. ; H. Cohl, 13 Dec, 1991 --- Added SOBEL keyword. ; H. Cohl, 7 Jun, 1991 --- Initial programming. ; T. Henry December, 1993 --- some changes made to size ;- function xyoff, a,b,bx,by,hx,hy,xx,yy,amax,fna,wdow=wdow,mask=mask,sobel=sobel, $ ccf=ccf,one=one sobel = keyword_set(sobel) ssz = size(a) ; Size of first image if ssz(0) ne 2 then message, 'Main image must be 2D array nx = ssz(1) ny = ssz(2) ;message, /info, 'Main image is A('+ $ ; strcompress(nx,/rem)+','+strcompress(ny,/rem)+')' ssz = size(b) if ssz(0) ne 2 then message, 'Secondary image must be 2D image if nx ne ssz(1) or ny ne ssz(2) then begin message, /info, 'Secondary image is A('+ $ strcompress(ssz(1),/rem)+','+strcompress(ssz(2),/rem)+')' message, 'Both images must have same array dimensions endif if n_elements(hx) eq 0 then hx = nx/2 ; Check center of image if n_elements(hy) eq 0 then hy = ny/2 if n_elements(bx) eq 0 then bx = nx-1 if n_elements(by) eq 0 then by = ny-1 mx = bx/2 & my = by/2 ; Set correlation box size. if not keyword_set(wdow) then wdow=1. ; Test to see if window is specified. if not keyword_set(mask) then mask=1. xlo = hx-mx & xhi = hx+mx ylo = hy-my & yhi = hy+my ; In the original XYOFF program the values of xhi and yhi were one less ; than in the current version. Activate the next line to go back to the ; original version ; xhi = xhi-1 & yhi = yhi-1 ;message, /info, 'Center of correlation window: ('+ $ ; strcompress(hx,/rem)+','+strcompress(hy,/rem)+')' ;message, /info, 'Window used for correlation: ' + $ ; 'A('+strcompress(xlo,/rem)+':'+strcompress(xhi,/rem)+','+ $ ; strcompress(ylo,/rem)+':'+strcompress(yhi,/rem)+')' if xlo lt 0 or xhi gt nx-1 or $ ylo lt 0 or yhi gt ny-1 then message, 'Invalid subarray index range na = a(xlo:xhi,ylo:yhi) if min(na) eq max(na) then message, 'First image is flat na = na-total(na)/n_elements(na) ; Subtract average if sobel then na = 1.*sobel(na) else na = na*wdow fna = fft(na*mask,-1) ; Fourier transform rms_na = sqrt(total(na*na)/n_elements(na)) nb = b(xlo:xhi,ylo:yhi) if min(na) eq max(na) then message, 'Second image is flat nb = nb-total(nb)/n_elements(nb) if sobel then nb = 1.*sobel(nb) else nb = nb*wdow fnb = fft(nb*mask,-1) ; Fourier transform rms_nb = sqrt(total(nb*nb)/n_elements(nb)) ccf = fft(fna*conj(fnb),1) ; Correlation fnc between na and nb ccf = float(ccf) ; Get real part of cross correlation fnc ccf = ccf/(rms_na*rms_nb) ; Normalize ccf = shift(ccf,mx,my) ; Shift correlation function. ps = max_pos(ccf) ; Find maximum of ccf s = [ps(0)-mx, ps(1)-my, max(ccf), ccf(bx/2,by/2), max(abs(ccf)) ] ;help, ccf ; xx = s(0) yy = s(1) amax = s(2); ; printf,one, xy,amax ; ;print, 'RMS IMAGE A & B : ', rms_na, rms_nb ;print, 'Position of Maximum(CCF)= ', [mx,my]+s(0:1) ;print, 'Shift in X - ', s(0) ;print, 'Shift in Y - ', s(1) ;print, 'Maximum(CCF) - ', s(2) ;print, 'Center CCF - ', s(3) ; ss=fltarr(9) ; ss(0) = ccf(bx/2-4,by/2) ; ss(1) = ccf(bx/2-3,by/2) ; ss(2) = ccf(bx/2-2,by/2) ; ss(3) = ccf(bx/2-1,by/2) ; ss(4) = ccf(bx/2,by/2) ; ss(5) = ccf(bx/2+1,by/2) ; ss(6) = ccf(bx/2+2,by/2) ; ss(7) = ccf(bx/2+3,by/2) ; ss(8) = ccf(bx/2+4,by/2) ;print,'x cor values ',ss(0),ss(1),ss(2),ss(3),ss(4),ss(5),ss(6),ss(7),ss(8) ;print, 'Maximum( ABS(CCF) )= ', s(4) ;print, 'Position of Maximum( ABS(CCF) )= ',max_pos(abs(ccf)) if n_elements(one) ne 0 then begin ;printf,one, 'RMS IMAGE A & B : ', rms_na, rms_nb ;printf,one, 'Position of Maximum(CCF)= ', [mx,my]+s(0:1) ;printf,one, 'Shift in X - ', s(0) ;printf,one, 'Shift in Y - ', s(1) ;printf,one, 'Maximum(CCF) - ', s(2) ;printf,one, 'Center CCF - ', s(3) ;printf,one,'x cor values ',ss(0),ss(1),ss(2),ss(3),ss(4),ss(5),ss(6),ss(7),ss(8) ;printf,one, 'Maximum( ABS(CCF) )= ', s(4) ;printf,one, 'Position of Maximum( ABS(CCF) )= ',max_pos(abs(ccf)) endif return, s & end ;------------------------------------------------------------------------- ;+ ;+ ; NAME: ; SFLOAT ; PURPOSE: ; Extract numbers from string and put them in floating point array ; CATEGORY: ; String processing ; CALLING SEQUENCE: ; FLTVEC = sfloat(INPUT,/exp,strvec=STRVEC,fmt=FMT, $ ; xfora=XFORA, numfmt=NUMFMT,lenfmt=LENFMT, $ ; strcrumbs=STRCRUMBS,numcrumbs=NUMCRUMBS,lencrumbs=LENCRUMBS) ; INPUTS: ; INPUT string (e.g. 89/1/12) ; EXP if set and nonzero, exponentials of type 'D' and 'E' are ; also interpreted ; OUTPUTS: ; FLTVEC floating array (e.g. [89,1,12]) ; STRVEC string array (array FLTVEC before conversion to floating point) ; FMT format string matching the input string ; NUMFMT # valid numbers located in the input string ; LENFMT # chars taken up by the NUMFMT numbers in the input string ; STRCRUMBS string array with part of input string not translated into ; numbers ; NUMCRUMBS # elements in STRCRUMBS ; LENCRUMBS # total length of all elements in STRCRUMBS summed ; PROCEDURE: ; Checks each character against a list of valid characters. Processing ; characters sequentially going from first to last, the largest possible ; substrings representing valid numbers are extracted. Only integer ; exponents are accepted (if the keyword EXP is set). ; MODIFICATION HISTORY: ; Written Jan 1989 by DMZ (ARC) ; Converted to V2 Dec 1990 by DMZ (ARC) ; OCT-1990, Paul Hick (ARC), Expanded to accept exponentials, ; ?-1993, Paul Hick (UCSD), Added construction of format specifiers ; FEB-1995, Paul Hick (UCSD), Added option to pull out the substrings ; not translated into numbers ;- ;----------------------------------------------------------------------- ;+ ; NAME: ; PCURSOR ; PURPOSE: ; Kludging the IDL CURSOR procedure ; CALLING SEQUENCE: ; PCURSOR, XP, YP, /device returns device coordinates ; PCURSOR, XP, YP, /normal returns normal coordinates ; PCURSOR, XP, YP, /data returns data coordinates ; FUNCTIONS/SUBROUTINES: ; CURSOR ; INPUTS: ; None ; OUTPUTS: ; XP, YP coordinates of spin-hair cursor ; !ERR = 1/2/4 when left/middle/right button is pushed ; PERROR !ERR value returned by CURSOR procedure ; KLUDGE if !D.NAME eq 'TEK' a blank print statement is inserted ; after the call to cursor ; COMMON BLOCKS: ; common DEVICES, PLOTDEV ; common SAVE_PCURSOR, GINCHARS ; KLUDGE: ; the KLUDGE keyword is required only for EMUTEK at the moment. Without ; the extra print statement the first call to CURSOR will work, but ; all subsequent calls will make EMUTEK get stuck in GIN mode. The ; reason is unclear. ; PROCEDURE: ; When using a mouse to select a screen location, the mouse buttons ; should return 1,2 and 4 in the !ERR system variable for left, middle ; and right mouse button (if the mouse has only two buttons use 1 and 2). ; The keys 1,2,3 or a,b,c or A,B,C or CTRL-A,CTRL-B,CTRL-C can be used ; to mimick mouse action, provided the CURSOR procedure puts the ASCII ; decimal equivalent of the corresponding key in the !ERR variable. If ; you get funny characters on the screen and/or you have to press one or ; more keys to return to the IDL prompt, there may be problems with the ; type-ahead buffer (i.e. CURSOR puts more/less characters in the type ; ahead buffer than IDL reads out of it before continuing). Set ; DEVICE,gin_chars=6 and try again. ; ; X-windows (!D.NAME = 'X'): ; Should work without any problems ; HDS3200 terminal (!D.NAME = 'TEK'): ; The mouse button keys should be user defined as , ; and for left, middle and right button respectively (the ; corresponding keyboard strokes are CTRL A, CTRL B and CTRL D). ; The error codes returned when pushing a button are the decimal ; equivalent of the first character of the user definition, which ; for the settings given above results 1,2 and 4 for left, middle ; and right button respectively. The IDL CURSOR procedure places ; the user definition of the button in the type-ahead buffer ; (one character in this case). ; VT240 (!D.NAME = 'REGIS'): ; It worked the last time I tried ; MODIFICATION HISTORY: ; FEB-91, Paul Hick (ARC) ;- pro PCURSOR, XP,YP, device=DEVICE,normal=NORMAL,data=DATA, $ perror=PERROR,kludge=KLUDGE,down=DOWN on_error, 1 ; On error return to main level common DEVICES, PLOTDEV common SAVE_PCURSOR, GINCHARS ;if !D.NAME eq PLOTDEV then return if !D.NAME eq 'TEK' then if n_elements(GINCHARS) eq 0 then begin GINCHARS = 6 & device, gin_chars=GINCHARS endif cursor, XP,YP,3,device=DEVICE, normal=NORMAL, data=DATA, down=DOWN PERROR = !ERR ; Get cursor location if !D.NAME ne 'X' then begin L = [1,49,65,97,76,108] ; Left button proxies: [SOH,1,A,a,L,l] M = [2,50,66,98,77,109] ; Middle button proxies: [STX,2,B,b,M,m] L = where(PERROR eq L) & M = where(PERROR eq M) if L(0) ne -1 then !ERR = 1 else if M(0) ne -1 then !ERR = 2 else !ERR = 4 endif if keyword_set(KLUDGE) and !D.NAME eq 'TEK' then print ; Kludge print statement for EMUTEK return & end ;----------------------------------------------------------------------- function SFLOAT, INPUT, exp=EXP, strvec=STRVEC, $ xfora=XFORA,fmt=FMT,numfmt=NUMFMT,lenfmt=LENFMT,$ strcrumbs=STRCRUMBS,numcrumbs=NUMCRUMBS,lencrumbs=LENCRUMBS if keyword_set(EXP) then EXP = 1 else EXP = 0 if keyword_set(XFORA) then XFORA = 1 else XFORA = 0 VALID_NR = ['1','2','3','4','5','6','7','8','9','0'] VALID_NRDOT = [VALID_NR,'.'] VALID_EXP = [VALID_NR,'+','-'] VALID_FULL = [VALID_NR,'.','+','-'] STRING = strtrim(INPUT) LEN = strlen(STRING) I = -1 & FLTVEC = '' & STRVEC = '' FMT = '*' NUMFMT = 0 FMTCNT = 0 LENFMT = 0 NUMCRUMBS = 0 LENCRUMBS = 0 LAST = I ; End of previous number BLANKS = 0 FIRST_CHAR: I = I+1 if I ge LEN then goto, DONE CH = strmid(STRING,I,1) if CH eq ' ' then begin BLANKS = BLANKS+1 goto, FIRST_CHAR endif A = where(CH eq VALID_FULL,CNUM) ; First char must be +,-,. or digit if CNUM eq 0 then begin BLANKS = 0 goto, FIRST_CHAR endif FIRST = -1 if CH eq '+' or CH eq '-' then begin J = I+1 ; Position following +,- if J ge LEN then begin I = J goto, DONE endif CH = strmid(STRING,J,1) A = where(CH eq VALID_NRDOT,CNUM) ; +,- must be followed by . or digit if CNUM eq 0 then begin BLANKS = 0 goto, FIRST_CHAR endif FIRST = I ; Position of +,- (start of number) I = J ; Position of . or digit endif if CH eq '.' then begin J = I+1 ; Position following . if J ge LEN then begin I = J goto, DONE endif CH = strmid(STRING,J,1) A = where(CH eq VALID_NR,CNUM) ; . must be followed by number if CNUM eq 0 then begin BLANKS = 0 goto, FIRST_CHAR endif DOT = I ; Position of . FMTCHAR = 'F' if FIRST eq -1 then FIRST = I ; Start of number I = J ; Position of digit endif else DOT = -1 if FIRST eq -1 then FIRST = I ; Start of number FIRST = FIRST-BLANKS ; Include blanks preceding number in format BLANKS = 0 ; Clear blanks counter if FIRST gt LAST+1 then begin case FMTCNT of 0: FMTL = 0 1: FMT = FMT+','+FMTSUB else: FMT = FMT+','+strcompress(FMTCNT,/remove_all)+FMTSUB endcase LENFMT = LENFMT+FMTCNT*FMTL A = FIRST-LAST-1 FMTCNT = 0 case XFORA of 0: FMT = FMT+',A'+strcompress(A,/remove_all) 1: FMT = FMT+','+strcompress(A,/remove_all)+'X' endcase NUMCRUMBS = NUMCRUMBS+1 LENCRUMBS = LENCRUMBS+A A = strmid(STRING,LAST+1,A) case NUMCRUMBS of 1: STRCRUMBS = A else: STRCRUMBS = [STRCRUMBS,A] endcase endif ; At this point the first digit is located (possibly preceded by +,+.,- or -.) ; i.e. a valid number has been located (as yet incomplete) ; FIRST = position where number starts ; DOT = position of . (-1 if no . was found) ; I = position of first digit ; Now start looking for the end of the number LAST_CHAR: I = I+1 if I ge LEN then goto, NUMBER CH = strmid(STRING,I,1) if DOT eq -1 and CH eq '.' then begin ; . is acceptable if it's the first one DOT = I FMTCHAR = 'F' goto, LAST_CHAR endif A = where(CH eq VALID_NR,CNUM) ; only digits are acceptable if CNUM gt 0 then goto, LAST_CHAR ; if the keyword EXP is not set then the number is considered completed. ; if the keyword EXP is set then continue only if the letter D or E is found ; If the number is complete then I is the position following the number if not EXP or (CH ne 'D' and CH ne 'E' and CH ne 'd' and CH ne 'e') then goto, NUMBER FMTCHAR = strupcase(CH) ; Analyze exponent IE = I ; Position of D,E IE = IE+1 if IE ge LEN then goto, NUMBER CH = strmid(STRING,IE,1) A = where(CH eq VALID_EXP,CNUM) ; D,E must be followed by +,- or number if CNUM eq 0 then goto, NUMBER if CH eq '+' or CH eq '-' then begin ; If IE is position of +,- J = IE+1 ; Position following +,- if J ge LEN then goto, NUMBER CH = strmid(STRING,J,1) A = where(CH eq VALID_NR,CNUM) ; +,- must be followed by number if CNUM eq 0 then goto, NUMBER IE = J ; Position of digit endif I = IE ; first digit in exponent located ; (possibly preceded by + or -) LAST_CHAR_EXP: I = I+1 if I ge LEN then goto, NUMBER CH = strmid(STRING,I,1) A = where(CH eq VALID_NR,CNUM) ; only numbers are acceptable if CNUM gt 0 then goto, LAST_CHAR_EXP NUMBER: STRVAL = strmid(STRING,FIRST,I-FIRST) FLVAL = float(STRVAL) if STRVEC(0) eq '' then begin STRVEC = STRVAL FLTVEC = FLVAL endif else begin STRVEC = [STRVEC,STRVAL] FLTVEC = [FLTVEC,FLVAL] endelse NUMFMT = NUMFMT+1 if DOT eq -1 then FMTCHAR = 'I' ;if DOT eq -1 then FMTCHAR = 'F' if DOT eq -1 then DOT = I-1 FMTNEW = FMTCHAR + strcompress(I-FIRST,/remove_all) case FMTCHAR of 'F': FMTNEW = FMTNEW + '.' + strcompress(I-DOT-1,/remove_all) 'D': FMTNEW = FMTNEW + '.' + strcompress(I-4-DOT-1,/remove_all) 'E': FMTNEW = FMTNEW + '.' + strcompress(I-4-DOT-1,/remove_all) else: FMTNEW = FMTNEW endcase if FMTCNT eq 0 then begin FMTSUB = FMTNEW FMTCNT = 1 FMTL = I-FIRST endif else if FMTNEW eq FMTSUB then $ FMTCNT = FMTCNT+1 $ else begin case FMTCNT of 1: FMT = FMT+','+FMTSUB else: FMT = FMT+','+strcompress(FMTCNT,/remove_all)+FMTSUB endcase LENFMT = LENFMT+FMTCNT*FMTL FMTSUB = FMTNEW FMTL = I-FIRST FMTCNT = 1 endelse I = I-1 LAST = I ; Pos of last character processed goto, FIRST_CHAR DONE: case FMTCNT of 0: FMTL = 0 1: FMT = FMT+','+FMTSUB else: FMT = FMT+','+strcompress(FMTCNT,/remove_all)+FMTSUB endcase LENFMT = LENFMT+FMTCNT*FMTL if I gt LAST+1 then begin A = I-LAST-1 case XFORA of 0: FMT = FMT+',A'+strcompress(A,/remove_all) 1: FMT = FMT+','+strcompress(A,/remove_all)+'X' endcase NUMCRUMBS = NUMCRUMBS+1 LENCRUMBS = LENCRUMBS+A A = strmid(STRING,LAST+1,A) case NUMCRUMBS of 1: STRCRUMBS = A else: STRCRUMBS = [STRCRUMBS,A] endcase endif FMT = strmid(FMT,2,999) FMT = '('+FMT+')' return, FLTVEC & end ;------------------------------------------------------------------------- ;+ ; NAME: ; SFLTARR ; PURPOSE: ; Read 2D float array from ASCII file ; CATEGORY: ; I/O ; CALLING SEQUENCE: ; Array = SFLTARR ( File, nx=NX, ny=NY, /exp, /stop, $ ; atleast=ATLEAST, scream=SCREAM, nocheck=NoCheck, $ ; timer=TIMER, fmts=FMTS, errormessage=errormessage, $ ; crumbs=Crumbs, comments=Comments ) ; FUNCTIONS/PROCEDURES: ; BOOST, SFLOAT ; INPUTS: ; FILE char name of ASCII file ; OPTIONAL INPUTS: ; nx=NX 1st dimension of output array Array (if omitted ; the number of elements in the 1st record is used) ; ny=NY 2nd dimension of output array Array (if omitted NY is ; set to the numbers of records in the file) ; atleast=ATLEAST only records with at least ATLEAST numbers in them are ; accepted ; /exp the records read from FILE are checked for exponents ; (see SFLOAT procedure) ; /stop only the first NY records are read ; /nocheck reduces double-checking in SfloatClean ; (see SfloatClean procedure) ; /skipfmt analyze each record with SFLOAT, i.e. do NOT try ; the accumulated formats used to read previous records ; /silent suppresses most of the informational messages ; /scream displays extra info about formats used to read file ; /timer displays timer information (testing purposes only) ; OUTPUTS: ; SFLTARR float 2D array of floating point numbers ; if something goes wrong SFLTARR=0 is returned ; Crumbs string 2D array containing the character substrings ; not converted into numbers ; errormessage = '' (null string) if the file is properly read ; = string describing error is something goes wrong ; OPTIONAL OUTPUT PARAMETERS: ; NX,NY dimensions of Array ; PROCEDURE: ; > Empty records and records starting with a semi-colon (;) are ignored ; > The function SFLOAT is used to convert a record read from FILE as ; a string into a floating point array ; > Only the first NX numbers of each record are stored in subsequent rows ; of Array. NX is explicitly given as input or, if omitted, is set equal ; to the number of elements in the first record of FILE ; > If less than NX numbers are found in a record, the corresponding ; row in Array is padded with zeros. ; > A buffer array of size NX by NY is initialized before reading ; the file, where NY is explicitly input or set to # records in the file. ; If FILE contains more than NY records, records following record NY are ; appended to Array by the BOOST procedure (this is slow for big arrays) ; unless the keyword STOP is set and non-zero. ; MODIFICATION HISTORY: ; OCT-1992, Paul Hick (UCSD) ; FEB-1995, Paul Hick (UCSD), added the 'crumbs' option function SFLTARR, File, nx=NX, ny=NY, exp=EXP, stop=STOP, $ atleast=ATLEAST, silent=SILENT, scream=SCREAM, $ NoCheck=NoCheck, timer=TIMER, fmts=FMTS, $ errormessage=errormessage, crumbs=Crumbs, $ skipfmt=SKIPFMT, comments=Comments ;- on_error, 2 ; On error, return to caller errormessage = '' if n_elements(ATLEAST) eq 0 then ATLEAST = 1 SCREAM = keyword_set(SCREAM) SILENT = keyword_set(SILENT) TIMER = keyword_set(TIMER) SKIPFMT= keyword_set(SKIPFMT) fi = file_search(File) if fi(0) eq '' then begin errormessage = 'File not found: '+File if not SILENT then message, /info, errormessage return, 0 endif message, /info , 'Processing file '+fi(0) if TIMER then stopwatch, /set ; Open file read only. Stop on error. openr, /get_lun, IU, fi(0), error=I if I ne 0 then begin errormessage = !ERR_STRING ; Open error: stop if not SILENT then message, /info, errormessage return, 0 endif on_ioerror, IOERROR ; Establish IO error handler ; Read header records until the first non-empty record which does not start ; with a semi-colon is found iREC = 0L iCOM = -1 STREC = '' ; STREC must exist for READF repeat begin if eof(IU) then message, 'Empty file '+fi(0); Empty file: stop readf, IU, STREC iREC = iREC+1 if strmid ( strtrim( STREC,1) ,0,1) eq ';' then begin if n_elements(Comments) eq 0 then $ Comments = STREC $ else $ Comments = [Comments,STREC] endif endrep until strcompress (STREC,/remove_all) ne '' and strmid ( strtrim( STREC,1) ,0,1) ne ';' ; Interpret 1st data record. Set array dimensions if not supplied as keyword. REC = sfloat(STREC,exp=EXP,fmt=FMT,numfmt=NUMFMT,lenfmt=LENFMT, $ strcrumbs=StrCrumbs, numcrumbs=nCrumbs) while NUMFMT lt ATLEAST do begin ; While no numbers in record if SCREAM then message, /info, $ 'Record'+strcompress(iREC)+' rejected: only'+strcompress(NUMFMT)+' numbers readf, IU, STREC iREC = iREC+1 REC = sfloat(STREC,exp=EXP,fmt=FMT,numfmt=NUMFMT,lenfmt=LENFMT, $ strcrumbs=StrCrumbs, numcrumbs=nCrumbs) endwhile NREC = n_elements(REC) if not keyword_set(NX) then NX = NREC ; X-dimension if not keyword_set(NY) then begin ; Y-dimension point_lun, -IU, FilePos ; Save file pointer NY = 1L ; # remaining records while not EOF(IU) do begin NY = NY+1 readf, IU, STREC endwhile point_lun, IU, FilePos ; Reset file pointer endif Array = fltarr(NX,NY) ; Create array if nCrumbs ne 0 then begin nCrumbsMax = nCrumbs Crumbs = strarr(nCrumbsMax,NY) endif else nCrumbsMax = 0 if not SILENT then message, /info, $ 'Created array A('+strcompress(NX,/rem)+','+strcompress(NY,/rem)+')' BOOSTED = 0 FMTS = FMT NUMSFMT = NUMFMT LENSFMT = LENFMT if SCREAM then message, /info, $ 'Record'+strcompress(iREC)+', new SFLOAT format : '+FMT ; Fill in 1st row in output array iROW = 0L & Ix = (n_elements(REC) < NX) - 1 ; Fill 1st row of array if iROW lt NY then begin Array(0:Ix,iROW) = REC(0:Ix) if nCrumbs*nCrumbsMax ne 0 then begin I = (nCrumbs < nCrumbsMax)-1 Crumbs(0:I,iROW) = StrCrumbs(0:I) endif endif else begin if not keyword_set(STOP) then begin boost, Array, REC(0:Ix) if nCrumbs*nCrumbsMax ne 0 then begin I = (nCrumbs < nCrumbsMax)-1 boost, Crumbs, StrCrumbs(0:I) endif NY = NY+1 BOOSTED = 1 endif else $ goto, DONE endelse on_ioerror, DIRECT_FAIL ; Establish IO error handler point_lun, -IU, FilePos ; Save file pointer ; Try using a single read to file the array (works hopefully only when ; all records match the same format FMT if SCREAM then message, /info, 'Reading file with format : '+FMT i = iROW+1 tmp = Array(*,i:*) readf, IU, form=FMT, tmp free_lun, IU Array(*,i:*) = tmp if NX lt NREC then message, /info, 'Output array is probably garbage. Try reading without NX keyword' if TIMER then stopwatch,/stop return, Array DIRECT_FAIL: on_ioerror, IOERROR ; Re-establish IO handler if not SILENT then message, /info, 'Error reading file in one big gulp: if not SILENT then message, /info, !ERR_STRING point_lun, IU, FilePos ; Reset file pointer ; Now comes the hard work: ; Process the remaining records until eof is reached while not eof(IU) do begin readf, IU, STREC iREC = iREC+1 if strcompress (STREC,/remove_all) ne '' and $ strmid ( strtrim(STREC,1) ,0,1) ne ';' then begin if SKIPFMT then goto, TRYSFLOAT on_ioerror, CONTINUE ; Loop over all accumulated formats for I=0,n_elements(FMTS)-1 do begin if strpos(FMTS(I),'A') ne -1 then $ sfloatclean, STREC, FMTS(I), STRECtmp, FMTtmp, NoCheck=NoCheck,$ strcrumbs=StrCrumbs, numcrumbs=nCrumbs $ else begin STRECtmp = STREC FMTtmp = FMTS(I) nCrumbs = 0 endelse if FMTtmp ne '()' then begin NUMFMT = NUMSFMT(I) REC = fltarr(NUMFMT,/nozero) reads, STRECtmp, format=FMTtmp, REC if SCREAM then message, /info, 'Record'+strcompress(iREC)+', read with : '+FMTS(I) on_ioerror, IOERROR goto, SKIP_SFLOAT ; If format fits, skip SFLOAT endif CONTINUE: endfor TRYSFLOAT: on_ioerror, IOERROR ; None of the formats fits; use SFLOAT REC = sfloat(STREC,exp=EXP,fmt=FMT,numfmt=NUMFMT,lenfmt=LENFMT, $ strcrumbs=StrCrumbs,numcrumbs=nCrumbs) I = where (FMT eq FMTS) if I(0) eq -1 then begin ; Save the format, if it's a new one FMTS = [FMTS, FMT] NUMSFMT = [NUMSFMT,NUMFMT] LENSFMT = [LENSFMT,LENFMT] I = rotate(sort(LENSFMT),5) ; Reverse the array: longest format first FMTS = FMTS(I) NUMSFMT = NUMSFMT(I) LENSFMT = LENSFMT(I) endif if SCREAM then message, /info, $ 'Record'+strcompress(iREC)+', new SFLOAT Format : '+FMT SKIP_SFLOAT: if NUMFMT ge ATLEAST then begin ; Record contains numbers REC = REC(0:(n_elements(REC) < NX)-1 ) iROW = iROW+1 & Ix = (n_elements(REC) < NX) - 1 if iROW lt NY then begin Array(0:Ix,iROW) = REC(0:Ix) if nCrumbs*nCrumbsMax ne 0 then begin I = (nCrumbs < nCrumbsMax)-1 Crumbs(0:I,iROW) = StrCrumbs(0:I) endif endif else begin if not keyword_set(STOP) then begin boost, Array, REC(0:Ix) if nCrumbs*nCrumbsMax ne 0 then begin I = (nCrumbs < nCrumbsMax)-1 boost, Crumbs, StrCrumbs(0:I) endif NY = NY+1 BOOSTED = 1 endif else $ goto, DONE endelse endif endif endwhile ; Close file, cancel I/O error handler ; Remove empty trailing rows in output array (if present) DONE: free_lun, IU & on_ioerror, NULL iROW = iROW+1 if iROW lt NY then begin Array = Array(*,0:iROW-1) if nCrumbsMax ne 0 then Crumbs = Crumbs(*,0:iROW-1) NY = iROW if not SILENT then message, /info, $ 'Array truncated to A('+strcompress(NX,/rem)+','+strcompress(NY,/rem)+')' endif else if BOOSTED then begin if not SILENT then message, /info, $ 'Array boosted to A('+strcompress(NX,/rem)+','+strcompress(NY,/rem)+')' endif if TIMER then stopwatch,/stop return, Array IOERROR: ; I/O error label errormessage = !ERR_STRING if not SILENT then message, /info, errormessage free_lun, IU & on_ioerror, NULL return, 0 & end ;-------------------------------------------------------------------------- ;+ ; NAME: ; RGBHelios ; PURPOSE: ; CATEGORY: ; CALLING SEQUENCE: ; FUNCTIONS/PROCEDURES: ; INPUTS: ; OPTIONAL INPUT PARAMETERS: ; OUTPUTS: ; OPTIONAL OUTPUT PARAMETERS: ; COMMON BLOCKS: ; SIDE EFFECTS: ; RESTRICTIONS: ; PROCEDURE: ; MODIFICATION HISTORY: ; DECLARATION: ;- pro MakeRGB on_error, 2 ;common devices , PlotDev,TermDev dname = !d.name ;set_plot, TermDev rr = 0 ; black gg = 0 bb = 0 ;loadct, 0 ; black-white: extract greyish ;tvlct, /get, r,g,b ;i = 60+4*indgen(40) ; 40 colors ;rr = [rr,r(i)] ;gg = [gg,g(i)] ;bb = [bb,b(i)] loadct, 1 ; blue-white: extract blueish tvlct, /get, r,g,b i = 100+3*indgen(40) ; 40 colors rr = [rr,r(i)] gg = [gg,g(i)] bb = [bb,b(i)] loadct, 5 ; std-gamma: extract yellow tvlct, /get, r,g,b i = [159+1*indgen(18),177+2*indgen(22)] ; 40 colors rr = [rr,r(i)] gg = [gg,g(i)] bb = [bb,b(i)] loadct, 3 ; red temperature: extract red tvlct, /get, r,g,b i = 80+3*indgen(40) ; 40 colors rr = [rr,r(i)] gg = [gg,g(i)] bb = [bb,b(i)] loadct, 8 ; green/white linear: extract green tvlct, /get, r,g,b i = 80+4*indgen(40) ; 40 colors rr = [rr,r(i)] gg = [gg,g(i)] bb = [bb,b(i)] i =[1,1] rr = [rr,252*i] ; purple (2 colors) gg = [gg, 76*i] bb = [bb,193*i] rr = [rr,255*i] ;red (2 colors) gg = [gg, 0*i] bb = [bb, 0*i] rr = [rr,102*i] ; 0*i] ;green (2 colors) gg = [gg,162*i] ; 255*i] bb = [bb, 0*i] ; 10*i] rr = [rr, 0*i] ;blue (2 colors) gg = [gg, 0*i] bb = [bb,255*i] ;map = read_bmp(filepath(root=getenv('sam'),'earth.bmp'),r,g,b) ;i = n_elements(rr) ;write_bmp, filepath(root=getenv('sam'),'earthtmp.bmp'),i+map,shift(r,i),shift(g,i),shift(b,i) rr = [rr, r(0:3)] ;4 colors gg = [gg, g(0:3)] bb = [bb, b(0:3)] rr = [rr,255] ;white (1 color) gg = [gg,255] bb = [bb,255] n = n_elements(rr) if n lt !d.n_colors then begin rr = [rr,replicate(255,!d.n_colors-n)] gg = [gg,replicate(255,!d.n_colors-n)] bb = [bb,replicate(255,!d.n_colors-n)] endif loadct, 4 ; blue/green/red/yellow tvlct, /get, r, g, b rr = [rr,r] ; append predefined color table gg = [gg,g] bb = [bb,b] rgb = [transpose(rr),transpose(gg),transpose(bb)] openw, /get_lun, iu, filepath(root=getenv('sam'),'rgbhelios.txt') printf, iu, rgb free_lun, iu set_plot, dname tvlct, rr(0:235), gg(0:235), bb(0:235) ; load the custom table return end ;--------------------------------------------------------------- pro SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish on_error, 2 ; Read color tables from file S = sfltarr(filepath(root=getenv('sam'),'rgbhelios.txt')) nS = n_elements(S)/(3*2) reds = S(0,nS:*) greens = S(1,nS:*) blues = S(2,nS:*) ; Set up rgbo vector ; Helios: OpaMin = 0.05, DenMin = 50 or 60 ; OpaMin = 0.10, DenMin = 30 or 40 (individual tiff8 file) ; IPS : OpaMin = 0.00, DenMin = 10 ; CME only: OpaMin=0.10, DenMin = 1. ; SPSM : OpaMin = 0.05, DenMin = 70 ; OpaMax = 0.40, DenMax = 130 ColMin = 130 ; Range of colors to be used ColMax = 225< (n_elements(reds)-1) reds = reds [ColMin:ColMax] greens = greens[ColMin:ColMax] blues = blues [ColMin:ColMax] ColMax = ColMax-ColMin ColMin = 0 OpaMin = .05 ; Opacity range to be used OpaMax = .3;.4 ;.5 OpaCap = .5 ; Range of density values to be mapped ; The minimum density is scaled to zero; ; The maximum density is scaled to RgbTop ; Only scaled densities in the range [DenMin,DenMax] are displayed. ; Color indices used for densities DenMin = 70; 10 DenMax = 130;200;130 DenCap = DenMax+1 DenTop = RgbTop DenVal = DenMin+indgen(DenMax-DenMin+1) x = (DenVal-DenMin)/float(DenMax-DenMin) DenCol = round( ColMin+(ColMax-ColMin)*x ) DenOpa = OpaMin+(OpaMax-OpaMin)*x DenOpa = DenOpa^2 rgbo = bytarr(256,4) rgbo(DenVal,0) = DenOpa*reds (DenCol) rgbo(DenVal,1) = DenOpa*greens(DenCol) rgbo(DenVal,2) = DenOpa*blues (DenCol) rgbo(DenVal,3) = DenOpa*255 if DenCap le DenTop then begin rgbo(DenCap:DenTop,0) = OpaCap*reds (ColMax) rgbo(DenCap:DenTop,1) = OpaCap*greens(ColMax) rgbo(DenCap:DenTop,2) = OpaCap*blues (ColMax) rgbo(DenCap:DenTop,3) = OpaCap*255 endif nS = nS-1 reds = S(0,0:nS) greens = S(1,0:nS) blues = S(2,0:nS) White = 255 Black = 0 ;Greyish = [ 1, 40] ;Blueish = [ 41, 80] ;Yellowish = [ 81,120] ;Reddish = [121,160] ;Greenish = [161,200] Blueish = [ 1, 40] Yellowish = [ 41, 80] Reddish = [ 81,120] Greenish = [121,160] Purple = 201 Red = 203 Green = 205 return & end ;------------------------------------------------------------- ;max_pos.pro: REMOVE THIS LINE ;------------------------------------------------------------- ;+ ; NAME: ; MAX_POS ; PURPOSE: ; This function, given an array with a single well ; defined maximum, finds the maximum position ; of that array to decimal accuracy using an interpolation ; technique to find the maximum of an array from, Niblack, W., ; "An Introduction to Digital Image Processing", p 139. ; CATEGORY: ; IMAGE PROCESSING ; CALLING SEQUENCE: ; max = max_pos(f) ; INPUTS: ; f = An array. ; type: array,any type ; KEYWORD PARAMETERS: ; OUTPUTS: ; max = An array containing the position of the maximum of an array. ; type: vector,floating point,fltarr(2) ; max(0) = The decimal position of the x-maximum. ; type: scalar,floating point ; max(1) = The decimal position of the y-maximum. ; type: scalar,floating point ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; H. Cohl, 20 Nov, 1990 --- Initial implementation. ;- ;------------------------------------------------------------- function max_pos,f,help=help ;Display idl header if help is required. if keyword_set(help) or n_params() lt 1 then begin get_idlhdr,'max_pos.pro' max=-1 goto,finishup endif ssz=size(f) fmax=max(f,loc) xmax=loc mod ssz(1) ymax=loc/ssz(1) ;A more complicated interpolation technique to find the maximum of an array. ;from, Niblack, W., "An Introduction to Digital Image Processing", p 139. if (xmax*ymax gt 0) and (xmax lt (ssz(1)-1)) and (ymax lt (ssz(2)-1)) then begin denom=fmax*2-f(xmax-1,ymax)-f(xmax+1,ymax) xfra=(xmax-.5)+(fmax-f(xmax-1,ymax))/denom denom=fmax*2-f(xmax,ymax-1)-f(xmax,ymax+1) yfra=(ymax-.5)+(fmax-f(xmax,ymax-1))/denom xmax=xfra ymax=yfra endif rmax=[xmax,ymax] finishup: return,rmax end ; function bound,ix,iy,xc,yc,scx,scy,az,ht ; ; Determine the azimuth and height from the input ix and iy ; radtd = 57.29578 ; scx2 = scx*scx scy2 = scy*scy az = atan((iy - yc),(ix - xc))*radtd ht = sqrt((iy - yc)*(iy - yc)*scy2 + (ix - xc)*(ix - xc)*scx2) if az lt 0.0 then az = 360.0 + az ; END ; function bestr,arr,amp,ang,nx,ny,xc,yc,scx,scy, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,br,brm1,brm2,naz minxx = 2 minyy = 2 ivv = miny ivp = maxy ihh = minx ihp = maxx ampr = fltarr(naz+1) amprm1 = fltarr(naz+1) amprm2 = fltarr(naz+1) amprm3 = fltarr(naz+1) ; brm1 = fltarr(naz+1) ; brm2 = fltarr(naz+1) ; brm3 = fltarr(naz+1) xr = fltarr(naz+1) xrm1 = fltarr(naz+1) xrm2 = fltarr(naz+1) xrm3 = fltarr(naz+1) yr = fltarr(naz+1) yrm1 = fltarr(naz+1) yrm2 = fltarr(naz+1) yrm3 = fltarr(naz+1) ampr(*) = 0.0 br(*) = 0.0 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin iaz = fix((azb - az2) + 0.5) if amp(ixx,iyy) gt ampr(iaz) then begin ampr(iaz)=amp(ixx,iyy) br(iaz) = arr(ixx,iyy) xr(iaz)=ixx yr(iaz)=iyy endif endif endif ENDFOR ENDFOR amprm1(*) = 0.0 brm1(*) = 0.0 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin iaz = fix((azb - az2) + 0.5) xdif = abs(ixx - xr(iaz)) ydif = abs(iyy - yr(iaz)) if amp(ixx,iyy) gt amprm1(iaz) $ and xdif gt minxx and ydif gt minyy then begin amprm1(iaz)=amp(ixx,iyy) brm1(iaz) = arr(ixx,iyy) xrm1(iaz)=ixx yrm1(iaz)=iyy endif endif endif ENDFOR ENDFOR amprm2(*) = 0.0 brm2(*) = 0.0 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin iaz = fix((azb - az2) + 0.5) xdif = abs(ixx - xr(iaz)) ydif = abs(iyy - yr(iaz)) xdif1 = abs(ixx - xrm1(iaz)) ydif1 = abs(iyy - yrm1(iaz)) xdif2 = abs(xr(iaz) - xrm1(iaz)) ydif2 = abs(yr(iaz) - yrm1(iaz)) if amp(ixx,iyy) gt amprm2(iaz) $ and xdif gt minxx and ydif gt minyy $ and xdif1 gt minxx and ydif1 gt minyy $ and xdif2 gt minxx and ydif2 gt minyy then begin amprm2(iaz)=amp(ixx,iyy) brm2(iaz) = arr(ixx,iyy) xrm2(iaz)=ixx yrm2(iaz)=iyy endif endif endif ENDFOR ENDFOR ; print, brm1, xrm1, yrm1 ; print, brm2, xrm2, yrm2 ; end ; function bestv,sp,amp,ang,nx,ny,xc,yc,scx,scy,xsize,ysize, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestsp0,bestsp1,bestsp2, $ bestamp0,bestamp1,bestamp2,bestang0,bestang1,bestang2 minxx = xsize/3 minyy = ysize/3 ivv = miny ivp = maxy ihh = minx ihp = maxx if ivv lt 0 then ivv = 0 if ivp ge ny then ivp = ny - 1 if ihh lt 0 then ihh = 0 if ihp ge nx then ihp = nx - 1 bx0 = -1 by0 = -1 ;print, az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestamp0,bestsp0,bestang0 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if amp(ixx,iyy) gt bestamp0 then begin bestamp0 = amp(ixx,iyy) bestsp0 = sp(ixx,iyy) bestang0 = ang(ixx,iyy) bx0 = ixx by0 = iyy endif endif endif ENDFOR ENDFOR ;print, az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestamp0,bestsp0,bestang0 if bx0 eq -1 then goto, endit bx1 = -1 by1 = -1 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin xdif = abs(ixx - bx0) ydif = abs(iyy - by0) if amp(ixx,iyy) gt bestamp1 $ and xdif gt minxx and ydif gt minyy then begin bestamp1 = amp(ixx,iyy) bestsp1 = sp(ixx,iyy) bestang1 = ang(ixx,iyy) bx1 = ixx by1 = iyy endif endif endif ENDFOR ENDFOR if bx1 eq -1 then goto, endit bx2 = -1 by2 = -1 FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if bx1 ne -1 and by1 ne -1 then begin xdif = abs(ixx - bx0) ydif = abs(iyy - by0) xdif1 = abs(ixx - bx1) ydif1 = abs(iyy - by1) xdif2 = abs(bx0 - bx1) ydif2 = abs(by0 - by1) if amp(ixx,iyy) gt bestamp2 $ and xdif gt minxx and ydif gt minyy $ and xdif1 gt minxx and ydif1 gt minyy $ and xdif2 gt minxx and ydif2 gt minyy then begin bestamp2 = amp(ixx,iyy) bestsp2 = sp(ixx,iyy) bestang2 = ang(ixx,iyy) bx2 = ixx by2 = iyy endif endif endif endif ENDFOR ENDFOR endit: ; end ; function cimage,i1,i2,nh1,nv1,nh2,nv2,xc1,yc1,xc2,yc2, $ scx1,scy1,scx2,scy2,pspeed,dtime,az1,az2,ht1,ht2, $ minx,maxx,miny,maxy,amaxs,i11s,i22s,ccfs,istar, $ pspeedxy,xcent,ycent,kk,itrigger,imgg,mggg, $ az11,az22,ht11,ht22,xsize,ysize,one, $ angl,sspeed,fspeed,smax,badspeed,badangle,badcorr, $ cspeed,cangl,cmax ; ; This function obtains sets of arrays of the proper size and displacement ; to be used in a correlation ; radtd = 57.29578 rs = 695980.0 ; Sun radius in km ; istep = 1 ; step increment between correlated data sets in x jstep = 1 ; step increment between correlated data sets in y xsize2p = 2*xsize+1 ; 2 times the x size of the array to correlate ysize2p = 2*ysize+1 ; 2 times the y size of the array to correlate isizeh = fix(xsize/2 + 1) ; half of maximum x value of correlation jsizeh = fix(ysize/2 + 1) ; half of maximum y value of correlation scxsq = scx1*scx1 scysq = scy1*scy1 ; azb1 = az1 htb1 = ht1 azb2 = az2 htb2 = ht2 azb11 = az11 htb11 = ht11 azb22 = az22 htb22 = ht22 ; ; ; Arrays into which to place correlated parameters ; cspeed = fltarr(nh1,nv1) cangl = fltarr(nh1,nv1) cmax = fltarr(nh1,nv1) ; ; Zero arrays of correlated values ; nh1m = nh1 - 1 nv1m = nv1 - 1 FOR i=0L,nh1m do BEGIN FOR j=0L,nv1m do BEGIN cspeed(i,j) = badspeed ; no correlation speed cangl(i,j) = badangle ; no angle between radial and corr. cmax(i,j) = badcorr ; no correlation value ENDFOR ENDFOR ; ; Correlation arrays ; i11 = fltarr(xsize2p,ysize2p) i22 = fltarr(xsize2p,ysize2p) is11 = fltarr(xsize2p,ysize2p) ; print, ' Enter correlation' ; ivvv = ysize if miny gt ivvv then ivvv = miny ; ivvp = nv1m - ysize if maxy lt ivvp then ivvp = maxy ; ihhh = xsize if minx gt ihhh then ihhh = minx ; ihhp = nh1m - xsize if maxx lt ihhp then ihhp = maxx ; ; Preview data within first boundary ; maxdif = 15.0 i1min = 100000.0 i1max = -100000.0 FOR iy=ivvv,ivvp,jstep do BEGIN FOR ix=ihhh,ihhp,istep do BEGIN ; ; Determine the boundary within which to do the correlation on image 1 ; bd1 = bound(ix,iy,xc1,yc1,scx1,scy1,azb,htb) if azb ge azb2 and azb le azb1 then begin if htb ge htb1 and htb le htb2 then begin ; ivv = iy-ysize ivp = iy+ysize ihh = ix-xsize ihp = ix+xsize ; ; Check differences in the first array to use in the correlation. ; iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN is11(ixh,iyv) = i1(ixx,iyy) ixh = ixh + 1 ENDFOR iyv = iyv + 1 ENDFOR diffmxmn = (max(is11(*,*)) - min(is11(*,*))) maxis11 = diffmxmn minis11 = diffmxmn if maxis11 gt i1max then i1max = maxis11 if minis11 lt i1min then i1min = minis11 endif endif ENDFOR ENDFOR if mggg eq 2 then begin st1dmaxd = i1min + maxdif endif else begin st1dmaxd= i1max endelse st2dmaxd = st1dmaxd print, 'i1max, i1min, st1dmaxd', i1max, i1min, st1dmaxd ; nxstar = 0L nxxx = 0L nxxxx = 0L kk = 0L FOR iy=ivvv,ivvp,jstep do BEGIN FOR ix=ihhh,ihhp,istep do BEGIN ; ; Determine the boundary within which to do the correlation on image 1 ; bd1 = bound(ix,iy,xc1,yc1,scx1,scy1,azbb,htbb) if azbb ge azb2 and azbb le azb1 then begin if htbb ge htb1 and htbb le htb2 then begin ; ivv = iy-ysize ivp = iy+ysize ihh = ix-xsize ihp = ix+xsize ; ; print, '1st ivv ivp ihh ihp', ivv,ivp,ihh,ihp ; ; Load up the first array to use in the correlation. ; iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN i11(ixh,iyv) = i1(ixx,iyy) ixh = ixh + 1 ENDFOR iyv = iyv + 1 ENDFOR ; iflag1 = 0 iflag2 = 0 iflag3 = 0 iflag4 = 1 ; ; Determine a new x, y position of the center of the second window to ; correlate by assuming radial outward motion at the preferred speed. ; phid = atan((iy - yc1)*scx1,(ix - xc1)*scy1) ; ; Preferred x displacement to next image ; ixd = fix((pspeed*dtime)/scx1*cos(phid)) ; ; Preferred y displacement to next image ; iyd = fix((pspeed*dtime)/scy1*sin(phid)) ; inx = ix + ixd ; inx = fix(inx) iny = iy + iyd ; iny = fix(iny) ; ; Determine if the boundary is within the boundary correlation of image 2 ; bd2 = bound(inx,iny,xc2,yc2,scx2,scy2,azb,htb) if azb ge azb22 and azb le azb11 then begin if htb ge htb11 and htb le htb22 then begin ; ; Load up the second array to use in the correlation. ; ivv = iny-ysize ivp = iny+ysize ihh = inx-xsize ihp = inx+xsize ; print, '2nd ivv ivp ihh ihp', ivv,ivp,ihh,ihp iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN if iyy lt 0 or iyy ge nv2 or ixx lt 0 or ixx ge nh2 $ then iflag4 = 0 else BEGIN i22(ixh,iyv) = i2(ixx,iyy) ; if i11(ixh,iyv) ne i22(ixh,iyv) then iflag1 = 1 if ixh ge isizeh and iyv ge jsizeh and $ ixh lt (xsize2p-isizeh) and iyv lt (ysize2p-jsizeh) then begin if i11(ixh,iyv) ne i11(isizeh,jsizeh) then iflag2 = 1 if i22(ixh,iyv) ne i22(isizeh,jsizeh) then iflag3 = 1 endif ; ixh = ixh + 1 endelse ENDFOR iyv = iyv + 1 ENDFOR ; endif endif ; ; The load-up of arrays to use in the correlation is now finished. ; Now, correlate the arrays as long as there is data in the arrays and ; as all the data in one of the arrays is not flat (without structure). ; iflag5 = 0 if iflag1 eq 1 then begin if iflag2 eq 1 then begin if iflag3 eq 1 then begin if iflag4 eq 1 then begin ; i1Xi2 = xyoff(i11,i22,xsize,ysize,xsize,ysize,xxi,yyj, $ amax,ccf=ccf,one=one) ; nxxx = nxxx + 1L ; ; print, 'xxi yyj amax', xxi,yyj,amax if abs(xxi) gt xsize/3 or abs(yyj) gt xsize/3 then begin xys = badspeed dang = badangle amax = badcorr endif else begin ; ; Now check to see if the images are starlike ; ; print, 'istar(0,0)', istar(0,0) ; mst = 0.6 ; If corr. better than this, the image is a star amax1s = 0.0 amax2s = 0.0 if istar(0,0) ne -10000.0 then begin ; ; Compare image 1 with a star ; st1dif = max(i11(*,*)) - min(i11(*,*)) ; print, 'st1dif st1dmaxd', st1dif, st1dmaxd if st1dif le st1dmaxd then begin ; if st1dif le 30.0 then begin i1Xst = xyoff(i11,istar,xsize,ysize,xsize,ysize,xxi,yyj, $ amax1s) endif else begin amax1s = 1.0 endelse ; ; Compare image 2 with a star ; st2dif = max(i22(*,*)) - min(i22(*,*)) ; print, 'st2dif st2dmaxd',st2dif, st2dmaxd if st2dif le st2dmaxd then begin ; if st2dif le 30.0 then begin i2Xst = xyoff(i22,istar,xsize,ysize,xsize,ysize,xxi,yyj, $ amax2s) endif else begin amax2s = 1.0 endelse endif ; if amax1s le mst and amax2s le mst then iflag5 = 1 ; image not a star ; if iflag5 eq 0 then nxstar = nxstar + 1 ; if iflag5 eq 0 then print, 'Image a star - amax1s amax2s mst', amax1s, amax2s, mst ; xx = xxi - ixd yy = yyj - iyd phi = atan((iy - yc1)*scy1,(ix - xc1)*scx1)*radtd iphflag = 0 if phi lt 0.0 then begin phi = phi + 360.0 iphflag = 1 endif the = atan(-yy,-xx)*radtd if the lt 0.0 then begin the = the + 360.0 if iphflag eq 0 then phi = phi + 360.0 endif dang = (the - phi) xy = sqrt(xx*xx*scxsq + yy*yy*scysq) xys = xy/dtime ; mini1 = min(i11(*,*)) ; maxi1 = max(i11(*,*)) ; print, 'mini1 maxi1 ix iy', mini1, maxi1, ix, iy ; mini2 = min(i22(*,*)) ; maxi2 = max(i22(*,*)) ; print,' mini2 maxi2 inx iny', mini2, maxi2, inx, iny ; print, 'xys dang amax',xys,dang,amax endelse ; print, yy,xx,xy,xys,iy,yc1,ix,xc1,phi*radtd,the*radtd ;if amax gt .9 then print, yy,xx,xy,xys,iy,yc,ix,xc,phi*radtd,the*radtd,'angle diff =',dang if abs(dang) gt angl then begin xys = badspeed dang = badangle amax = badcorr endif if amax lt smax then begin xys = badspeed dang = badangle amax = badcorr endif if xys lt sspeed then begin xys = badspeed dang = badangle amax = badcorr endif if xys gt fspeed then begin xys = badspeed dang = badangle amax = badcorr endif endif endif endif endif if iflag1 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag2 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag3 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag4 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag5 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif ; if xys ne badspeed and dang ne badangle and amax ne badcorr then begin nxxxx = nxxxx + 1L if amax gt amaxs then begin amaxs = amax i11s(*,*) = i11(*,*) i22s(*,*) = i22(*,*) ccfs(*,*) = ccf(*,*) print, 'Max amp image - amax1s amax2s amaxs xys dang azbb htbb', $ amax1s, amax2s, amaxs, xys, dang, azbb, htbb/rs endif ; endif cspeed(ix,iy) = xys cangl(ix,iy) = dang cmax(ix,iy) = amax ; ; if iflag1 eq 1 and iflag2 eq 1 and iflag3 eq 1 then begin if xys ne badspeed and dang ne badangle and amax ne badcorr then begin pspeedxy(kk) = cspeed(ix,iy) xcent(kk) = ix ycent(kk) = iy kk = kk + 1L endif ; endif endif ENDFOR ENDFOR ; if nxxx gt 0L then begin print, ' success ', nxxx, ' pixels ',xsize,' by ',ysize,' were correlated' printf,one,' success ', nxxx, ' pixels ',xsize,' by ',ysize,' were correlated' print, ' ', nxstar, ' pixels were stellar correlations' printf, one, ' ', nxstar, ' pixels were stellar correlations' if nxxxx gt 0L then begin print, ' more success ', nxxxx, ' pixels past all criterion' printf, one, ' more success ', nxxxx, ' pixels past all criterion' endif endif if nxxx eq 0L then begin print, ' failure - no pixels correlated printf, one,' failure - no pixels correlated endif END ; ; function findmx, az1,az2,ht1,ht2,xc,yc,scx,scy,xsize,ysize,iave, $ minx,maxx,miny,maxy ; ; This function finds the minimum and maximum x and y values from the values ; of azimuth and height. ; radtd = 57.29578 ; naz = abs(az1-az2) nazp = naz+1 nazp2 = nazp*2 htx = fltarr(nazp2) hty = fltarr(nazp2) htx(*) = -999 hty(*) = -999 ; for i=0L,nazp do begin si = sin((az1-i)/radtd) co = cos((az1-i)/radtd) htx(i) = xc + (ht1/scx)*co hty(i) = yc + (ht1/scy)*si htx(nazp2-1-i) = xc + (ht2/scx)*co hty(nazp2-1-i) = yc + (ht2/scy)*si endfor iav = xsize if ysize gt iav then iav = ysize if iave gt iav then iav = iave minx = min(htx) - iav minx = fix(minx) maxx = max(htx) + iav + 1.0 maxx = fix(maxx) miny = min(hty) - iav miny = fix(miny) maxy = max(hty) + iav + 1.0 maxy = fix(maxy) ; END ; function filtr,ar,nh,nv,iave,xsize,ysize,xc,yc,scx,scy, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,i1 ; ; This function filters an image file a within the boundaries of the windowed ; area ; radtd = 57.29578 ; i1 = fltarr(nh,nv) ; FOR iyy=0L,nv-1 do BEGIN FOR ixx=0L,nh-1 do BEGIN i1(ixx,iyy) = -5.0 ENDFOR ENDFOR ; nxxx = 0L ; iavx = iave if xsize gt iave then iavx = xsize iavy = iave if ysize gt iave then iavy = ysize iavs = iavx*scx if iavy*scy gt iavs then iavs = iavy*scy ; dazbd = abs(atan(iavs,ht1))*radtd ; azb1 = az1 + dazbd htb1 = ht1 - iavs azb2 = az2 - dazbd htb2 = ht2 + iavs ; print, 'iavs azb1 azb2 htb1 htb2', iavs,azb1,azb2,htb1,htb2 ; ivv = iave-1 if miny gt ivv then ivv = miny ; ivp = nv - iave if maxy lt ivp then ivp = maxy ; ihh = iave-1 if minx gt ihh then ihh = minx ; ihp = nh - iave if maxx lt ihp then ihp = maxx ; print, ' Enter filtering' FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) ; print, 'ixx iyy azb htb', ixx,iyy,azb,htb if azb ge azb2 and azb le azb1 then begin if htb ge htb1 and htb le htb2 then begin nxxx = nxxx + 1L xyave = ar(ixx,iyy) xxyy = sqrt((iyy -yc)*(iyy -yc) + (ixx -xc)*(ixx -xc)) if xxyy ne 0.0 then BEGIN ax = (ixx - xc)/xxyy endif else BEGIN ax = 0.0 endelse if xxyy ne 0.0 then BEGIN ay = (iyy - yc)/xxyy endif else BEGIN ay = 0.0 endelse FOR i=0L,iave-1 do BEGIN iix = i*ax iiy = i*ay ixm = ixx - iix iym = iyy - iiy ixp = ixx + iix iyp = iyy + iiy xyave = xyave + ar(ixm,iym) + ar(ixp,iyp) ENDFOR xyave = xyave/(2.0*iave-1) i1(ixx,iyy) = ar(ixx,iyy) - xyave endif endif ENDFOR ENDFOR print, ' success - filtering ',nxxx,' points' end ; function mxmnbv, a,xm,ym,az1,az2,ht1,ht2,xc,yc,scx,scy, $ minx,maxx,miny,maxy,mxbv,mnbv ; ; This function finds the minimum and maximum image values within the area ; bound by azimuth and height. ; radtd = 57.29578 ; mxbv = -100000.0 mnbv = 100000.0 for i=minx,maxx do begin for j=miny,maxy do begin if i gt 0 and i lt xm then begin if j gt 0 and j lt ym then begin bd = bound(i,j,xc,yc,scx,scy,azb,htb) if azb ge az2 and azb le az1 then begin if htb ge ht1 and htb le ht2 then begin if a(i,j) gt mxbv then mxbv = a(i,j) if a(i,j) lt mxbv then mnbv = a(i,j) endif endif endif endif endfor endfor END ; function smoo,ar,nh,nv,xsize,ysize,xc,yc,scx,scy, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,i1 ; ; This function smooths an image file within the boundaries of the windowed ; area so that point glitches, cosmic rays, especially are removed from ; the images. ; radtd = 57.29578 ; i2 = fltarr(nh,nv) i1 = fltarr(nh,nv) ; i1(*,*) = -5.0 i2(*,*) = ar(*,*) ; difmax = 800. ; iavx = xsize iavy = ysize ; iavs = iavx*scx if iavy*scy gt iavs then iavs = iavy*scy ; dazbd = abs(atan(iavs,ht1))*radtd ; azb1 = az1 + dazbd htb1 = ht1 - iavs azb2 = az2 - dazbd htb2 = ht2 + iavs ; ivv = iavy-1 if miny gt ivv then ivv = miny ; ivp = nv - iavy if maxy lt ivp then ivp = maxy ; ihh = iavx-1 if minx gt ihh then ihh = minx ; ihp = nh - iavx if maxx lt ihp then ihp = maxx ; print, ' Enter Smoothing' ; ; First pass ; nx1 = 0L FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) if azb ge azb2 and azb le azb1 then begin if htb ge htb1 and htb le htb2 then begin arsur = (ar(ixx+1)*2.0 + ar(ixx-1)*2.0 + ar(iyy+1)*2.0 $ + ar(iyy-1)*2.0 + ar(ixx+1,iyy+1) + ar(ixx-1,iyy-1) $ + ar(ixx+1,iyy-1) + ar(ixx-1,iyy+1))/12.0 ardiff = abs(ar(ixx,iyy) - arsur) if ardiff gt difmax then begin i2(ixx,iyy) = arsur nx1 = nx1 + 1L endif endif endif ENDFOR ENDFOR if nx1 gt 0 then print, 'First pass found', nx1, ' points to smooth' ar(*,*) = i2(*,*) ; ; Second pass ; nx2 = 0L FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) if azb gt azb2 and azb lt azb1 then begin if htb gt htb1 and htb lt htb2 then begin arsur = (ar(ixx+1)*2.0 + ar(ixx-1)*2.0 + ar(iyy+1)*2.0 $ + ar(iyy-1)*2.0 + ar(ixx+1,iyy+1) + ar(ixx-1,iyy-1) $ + ar(ixx+1,iyy-1) + ar(ixx-1,iyy+1))/12.0 ardiff = abs(ar(ixx,iyy) - arsur) if ardiff gt difmax then begin i2(ixx,iyy) = arsur nx2 = nx2 + 1L endif endif endif ENDFOR ENDFOR if nx2 gt 0 then print, 'Second pass found', nx2, ' points to smooth' ; i1(*,*) = i2(*,*) ; end ; function getcp,xc,yc,scx,scy,az1,az2,ht1,ht2,iok ; ; Determine an image area to analyze in the study. ; radtd = 57.29578 ; if iok eq 0 then begin ; pcursor, xp11, yp11, /DEVICE pcursor, xp22, yp22, /DEVICE ; print, 'xp11 yp11 xp22 yp22', xp11,yp11,xp22,yp22 xp1 = xp11 yp1 = yp11 xp2 = xp22 yp2 = yp22 scx2 = scx*scx scy2 = scy*scy az11 = atan((yp1 - yc),(xp1 - xc))*radtd if az11 lt 0.0 then az11 = az11 + 360.0 az22 = atan((yp2 - yc),(xp2 - xc))*radtd if az22 lt 0.0 then az22 = az22 + 360.0 print, az11, az22 yp1m = yp1 - yc yp2m = yp2 - yc if az11 lt az22 then az11 = 360.0 + az11 ht11 = sqrt((yp1 - yc)*(yp1 - yc)*scy2 + (xp1 - xc)*(xp1 - xc)*scx2) ht22 = sqrt((yp2 - yc)*(yp2 - yc)*scy2 + (xp2 - xc)*(xp2 - xc)*scx2) az1 = az11 az2 = az22 ht1 = ht11 ht2 = ht22 if ht1 gt ht2 then begin ht2 = ht11 ht1 = ht22 endif endif ; END ; function getok, ih, iv, okok ; ; If the cursor is clicked within the dot on the overplot box of imagea, then ; a "true" is returned for the function. ; getok = 1 ; the only time getok is false is if ; you say so with the cursor if okok eq 0 then begin pcursor, xp, yp, /DEVICE getok = xp gt ih and xp lt ih+20 and yp gt 0 and yp lt 20 if getok eq 1 then print, ' Keep this one!' if getok eq 0 then print, ' Try again!' endif return, getok END ; function getspeed,xc,yc,scx,scy,ht11,ht22,dt,psp,iok ; ; Determine an image area to analyze in the study. ; radtd = 57.29578 ; if iok eq 0 then begin ; ; Ht from cursor mark ; print, 'define speed - click on center of new box - speed is now',psp pcursor, xp11, yp11, /DEVICE ; xp1 = xp11 yp1 = yp11 scx2 = scx*scx scy2 = scy*scy ht2 = sqrt((yp1 - yc)*(yp1 - yc)*scy2 + (xp1 - xc)*(xp1 - xc)*scx2) ; ; Ht of earlier box center ; ht1 = (ht11 + ht22)/2.0 ; ; Preferred speed ; psp = (ht2 - ht1)/dt ; endif ; END ; function image,name,ih,iv,mxbv,mnbv,imgg,a ; ; Image reader ; this function reads in a LASCO difference image and plots it ; a = readfits(name,hdr) ; mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 print, 'mxbv mnbv mx mn', mxbv, mnbv, mx, mn if imgg eq 1 then begin MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = mx, $ min = mn, top = Reddish(1)-Reddish(0)) ; print, reddish(0),reddish(1),greenish(0),greenish(1),blueish(0),blueish(1) WINDOW, /free, xsize=ih+20, ysize=iv plot, [0,ih+19], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih+19], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 TV, a0 endif ; END ; function image3,spv,ih,iv,sspeed,fspeed,minc,maxc ; ; This function images the speed two-dimensional image inside the boundary. ; All other values are set to the slowest speed. ; WINDOW, /free, xsize=ih, ysize=iv plot, [0,ih-1], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih-1], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 a0 = minc + bytscl( spv, max = fspeed, $ min = sspeed, top = maxc-minc) TV, a0 ; END ; function image4,ang,ih,iv,angl,sspeed,fspeed,smax,az1,az2,ht1,ht2 ; ; This function images the angle two-dimensional image ; MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Greenish(0) + bytscl( a, max = angl, $ min = -angl, top = Greenish(1)-Greenish(0)) WINDOW, /free, xsize=ih, ysize=iv plot, [0,ih-1], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih-1], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 TV, a0 ; END ; function image5,sma,ih,iv,angl,sspeed,fspeed,smax,az1,az2,ht1,ht2 ; ; This function images the correlation maximum two-dimensional image ; MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = 1.0, $ min = 0.0, top = Reddish(1)-Reddish(0)) WINDOW, /free, xsize=ih, ysize=iv plot, [0,ih-1], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih-1], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 TV, a0 ; END ; function image9,name,spv,ang,sma,ih,iv,badspv,badang,badsma,one,imgg, $ xsiz,ysiz,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,mxbv,mnbv ; ; This function reads in a LASCO difference image, plots it and overplots ; arrows which map the outward speed of features within the window. ; if imgg eq 1 then begin radtd = 57.29578 ; xsz = xsiz/3 ; nearest x corr. step to arrow base ysz = ysiz/3 ; nearest y corr. step to arrow base ihm = ih - 1 ivm = iv - 1 ; a = readfits(name,hdr) ; WINDOW, /free, xsize = ih, ysize = iv plot, [0,ih-1], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih-1], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 print, 'mxbv mnbv mx mn', mxbv, mnbv, mx, mn MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; a0 = Reddish(0) + bytscl( a, max = mx, $ min = mn, top = Reddish(1)-Reddish(0)) WINDOW, /free, xsize=ih+20, ysize=iv plot, [0,ih+19], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih+19], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 TV, a0 ; pix = 10.0/100.0 ; arrow length scale (pixels per km/s) ahead = 5.0 ; length of arrowhead dphi = 25.0/radtd ; angular shift for arrowhead spvs = 0.0 angs = 0.0 smas = 0.0 nsums = 0.0 arx = fltarr(2) ary = fltarr(2) arhx = fltarr(3) arhy = fltarr(3) FOR i=0L, ihm, xsz DO BEGIN FOR j=0L, ivm, ysz DO BEGIN if spv(i,j) ne badspv and ang(i,j) ne badang and sma(i,j) ne badsma $ then begin print, ' success', i,j,spv(i,j),ang(i,j),sma(i,j) printf, one,' success', i,j,spv(i,j),ang(i,j),sma(i,j) spvs = spvs + spv(i,j) angs = angs + ang(i,j) smas = smas + sma(i,j) nsums = nsums + 1.0 phid = atan((j - vcent1)*scy1,(i - hcent1)*scx1)*radtd if phid lt 0.0 then phid = phid + 360.0 phi = (phid + ang(i,j))/radtd ; phi in radians xd = spv(i,j)*pix*cos(phi) yd = spv(i,j)*pix*sin(phi) xp = i + xd yp = j + yd xpm = xp - ahead*cos(phi-dphi) ypm = yp - ahead*sin(phi-dphi) xpp = xp - ahead*cos(phi+dphi) ypp = yp - ahead*sin(phi+dphi) ; ; set up plotting arrays ; arx(0) = i ary(0) = j arx(1) = xp ary(1) = yp arhx(0) = xpm arhy(0) = ypm arhx(1) = xp arhy(1) = yp arhx(2) = xpp arhy(2) = ypp ; oplot, arx, ary, thick=2, /noclip oplot, arhx, arhy, thick=2, /noclip endif ENDFOR ENDFOR if nsums ne 0.0 then begin spvs = spvs/nsums angs = angs/nsums smas = smas/nsums print, ' summary of successes ', nsums,spvs,angs,smas printf, one, ' summary of successes ', nsums,spvs,angs,smas print, ' ' print, ' ' printf, one, ' ' printf, one, ' ' endif endif ; END ; function image10,name,ih,iv,az1,az2,ht1,ht2,xc,yc,scx,scy,imgg, $ minx,maxx,miny,maxy,value,minva,maxva,mxbv,mnbv,minc,maxc ; ; Image reader ; ; This function reads in a LASCO difference image and plots it outside the ; boundary defined before. Inside the boundary is plotted the appropriate ; parameter - pixel by pixel. ; if imgg eq 1 then begin a = readfits(name,hdr) ; azb1 = az1 htb1 = ht1 azb2 = az2 htb2 = ht2 ; ivv = miny ivp = maxy ihh = minx ihp = maxx ; MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 a0 = Reddish(0) + bytscl( a, max = mx, $ min = mn, top = Reddish(1)-Reddish(0)) WINDOW, /free, xsize=ih, ysize=iv plot, [0,ih-1], [0,iv-1], /nodata, /data, XMARGIN = [0,0], YMARGIN =[0,0], $ XRANGE = [0,ih-1], YRANGE = [0,iv-1], XSTYLE = 5, YSTYLE = 5 TV, a0 dif = maxc - minc FOR iyy=ivv,ivp do BEGIN FOR ixx=ihh,ihp do BEGIN bd = bound(ixx,iyy,xc,yc,scx,scy,azb,htb) if azb gt azb2 and azb lt azb1 then begin if htb gt htb1 and htb lt htb2 then begin va = value(ixx,iyy) if dif gt 39 then begin if va lt 0 then begin vlu = minc + bytscl( [va], max = 0.0, min = minva, $ top = 39) endif if va ge 0 then begin vlu = maxc - bytscl( [va], max = maxva, min = 0.0, $ top = 39) endif endif else begin vlu = minc + bytscl( [va], max = maxva, min = minva, $ top = maxc-minc) endelse TV, vlu, ixx, iyy endif endif ENDFOR ENDFOR ; endif ; END ; function imagea,ih,iv,xc,yc,scx,scy,az1,az2,ht1,ht2,imgg, $ xsiz,ysiz,az1i,ht2i,ti1,ti2,dtime,ps ; ; this overplots the "area to be analyzed lines" on the image as well as an ; "OK" box on the bottom right of the image. ; if imgg eq 1 then begin ; radtd = 57.29578 ; naz = abs(az1-az2) htx1 = fltarr(naz+1) hty1 = fltarr(naz+1) htx2 = fltarr(naz+1) hty2 = fltarr(naz+1) ; for i=0L,naz do begin si = sin((az1-i)/radtd) co = cos((az1-i)/radtd) htx1(i) = xc + (ht1/scx)*co hty1(i) = yc + (ht1/scy)*si htx2(i) = xc + (ht2/scx)*co hty2(i) = yc + (ht2/scy)*si endfor ; azxp1 = fltarr(2) azyp1 = fltarr(2) azxp2 = fltarr(2) azyp2 = fltarr(2) ; abox = fltarr(5) aboy = fltarr(5) ; aboxs = fltarr(5) aboys = fltarr(5) ardisx = fltarr(2) ardisy = fltarr(2) ; azxp1(0) = htx1(0) azyp1(0) = hty1(0) azxp1(1) = htx2(0) azyp1(1) = hty2(0) azxp2(0) = htx1(naz) azyp2(0) = hty1(naz) azxp2(1) = htx2(naz) azyp2(1) = hty2(naz) ; abox(0) = ih+2 abox(1) = ih+2 abox(2) = ih+18 abox(3) = ih+18 abox(4) = ih+2 aboy(0) = 2 aboy(1) = 18 aboy(2) = 18 aboy(3) = 2 aboy(4) = 2 ; dsp = ps*(ti2-ti1) ; preferred distance travelled ; dsp100 = 100.*(dtime) ; 100 km/s distance travelled pdistx = xc + (ht2i+dsp)/scx*cos(az1i/radtd) ; preferred x distance pdisty = yc + (ht2i+dsp)/scy*sin(az1i/radtd) ; preferred y distance aboxs(0) = pdistx aboxs(1) = pdistx xsi = xsiz if az1i gt 90. then xsi = -xsiz aboxs(2) = pdistx + xsi aboxs(3) = pdistx + xsi aboxs(4) = pdistx aboys(0) = pdisty ysi = ysiz if az1i gt 180. then ysi = -ysi aboys(1) = pdisty + ysi aboys(2) = pdisty + ysi aboys(3) = pdisty aboys(4) = pdisty ; ardisx(0) = pdistx + xsi ; ardisx(1) = pdistx + xsi + dsp100/scx*cos(az1i/radtd) ; ardisy(0) = pdisty + ysi ; ardisy(1) = pdisty + ysi + dsp100/scy*sin(az1i/radtd) ; oplot, azxp1, azyp1, thick = 2, /noclip oplot, azxp2, azyp2, thick = 2, /noclip oplot, htx1, hty1, thick = 2, /noclip oplot, htx2, hty2, thick = 2, /noclip oplot, abox,aboy, thick = 2, /noclip oplot, aboxs,aboys, thick = 2, /noclip ; oplot, ardisx,ardisy, thick = 2, /noclip ; endif ; END ; function imageb,ih,iv,xc,yc,scx,scy,az1,az2,ht1,ht2,imgg, $ xsiz,ysiz,az1i,ht2i,ti1,ti2,dtime,ps, $ minv,maxv,minc,maxc ; ; this overplots the "area to be analyzed lines" on the image as well as a ; color label box at the bottom right of the image. ; if imgg eq 1 then begin radtd = 57.29578 ; naz = abs(az1-az2) htx1 = fltarr(naz+1) hty1 = fltarr(naz+1) htx2 = fltarr(naz+1) hty2 = fltarr(naz+1) ; for i=0L,naz do begin si = sin((az1-i)/radtd) co = cos((az1-i)/radtd) htx1(i) = xc + (ht1/scx)*co hty1(i) = yc + (ht1/scy)*si htx2(i) = xc + (ht2/scx)*co hty2(i) = yc + (ht2/scy)*si endfor ; azxp1 = fltarr(2) azyp1 = fltarr(2) azxp2 = fltarr(2) azyp2 = fltarr(2) ; abox = fltarr(5) aboy = fltarr(5) ; aboxs = fltarr(5) aboys = fltarr(5) ardisx = fltarr(2) ardisy = fltarr(2) ; azxp1(0) = htx1(0) azyp1(0) = hty1(0) azxp1(1) = htx2(0) azyp1(1) = hty2(0) azxp2(0) = htx1(naz) azyp2(0) = hty1(naz) azxp2(1) = htx2(naz) azyp2(1) = hty2(naz) ; ht = maxc - minc abox(0) = ih-18 abox(1) = ih-18 abox(2) = ih-2 abox(3) = ih-2 abox(4) = ih-18 aboy(0) = 20 aboy(1) = ht + 20 aboy(2) = ht + 20 aboy(3) = 20 aboy(4) = 20 ; scale1 = fltarr(16,ht) FOR j=0L,ht-1 do BEGIN FOR i=0L,15 do BEGIN scale1(i,j) = minc + j if j gt 39 then scale1(i,j) = maxc - j + 39 endfor ; print, ht, j, minc + j, maxc - j + 39 endfor ; dsp = ps*(ti2-ti1) ; preferred distance travelled pdistx = xc + (ht2i+dsp)/scx*cos(az1i/radtd) ; preferred x distance pdisty = yc + (ht2i+dsp)/scy*sin(az1i/radtd) ; preferred y distance ; aboxs(0) = pdistx aboxs(1) = pdistx xsi = xsiz if az1i gt 90. then xsi = -xsiz aboxs(2) = pdistx + xsi aboxs(3) = pdistx + xsi aboxs(4) = pdistx aboys(0) = pdisty ysi = ysiz if az1i gt 180. then ysi = -ysi aboys(1) = pdisty + ysi aboys(2) = pdisty + ysi aboys(3) = pdisty aboys(4) = pdisty ; scl = minc + bytscl( scale1, max = maxc, $ min = minc, top = maxc-minc) TV, scl, ih-18,20 oplot, azxp1, azyp1, thick = 2, /noclip oplot, azxp2, azyp2, thick = 2, /noclip oplot, htx1, hty1, thick = 2, /noclip oplot, htx2, hty2, thick = 2, /noclip oplot, abox,aboy, thick = 2, /noclip oplot, aboxs,aboys, thick = 2, /noclip ; endif ; END function imagec,hsize,vsize,minv,maxv,minc,maxc,imgg ; ; This function plots the scale on the box at the left of the image. ; if imgg eq 1 then begin ; if minv ge 0.0 then min = fix(minv) if minv lt 0.0 then min = - fix(abs(minv)) max = fix(maxv) dif = maxc - minc x1 = hsize - 85 y1 = 17 y2 = dif + 17 xyouts, x1, y1, min xyouts, x1, y2, max ; endif ; END ; function nextim,az1,az2,ht1,ht2,filta,filtb1,filtb2, $ azd,hts,htf,az11,az22,ht11,ht22 ; ; This function determines the "area to be analyzed locations" on the second ; image ; az11 = az1 + azd + filta az22 = az2 + azd - filta ht11 = ht1 + hts - filtb1 ; ht in km, time in s. ht22 = ht2 + htf + filtb2 ; ht in km, time in s. ; end ; function param,name,img,istep,fini1,name1,h1,v1,hc,vc,scx,scy,doy,tim ; ; This function reads the IDL data file and for this image, determines ; its important parameters and returns the image. ; ; name of differenced image ; ; ir = lasco_readfits(name1,hdr) ; namep = namep+strcompress(image,/rem)+finif ; namep = name image = img+istep name1 = namep+strcompress(image,/rem)+fini1 a = lasco_readfits(name1,hdr) ; h1 = hdr.NAXIS1 ; NAXIS1 horizontal size of differenced image v1 = hdr.NAXIS2 ; NAXIS2 vertical size of differenced image sun = getsuncen(hdr) hc = sun.xcen vc = sun.ycen scxas = hdr.PLATESCL ; CDELT1 horizontal scale of image (arcsec/pixel) scyas = hdr.PLATESCL ; CDELT2 vertical scale of image (arcsec/pixel) doy = hdr.MID_DATE ; MID_DATE time of image in Modified Julian Day tim = hdr.MID_TIME ; MID_TIME time of day of image in miliseconds ; ; Determine rsun = sun size in arc seconds ; date = {mjd:hdr.MID_DATE,time:hdr.MID_TIME*1000.} rsu = sun_ephem(date) rsun = rsu(4)*60.0 ; convert to size in arcseconds print, rsun ; rs = 695980.0 ; Sun radius in km kmpas = rs/rsun ; Km on sun per arcsec scx = kmpas*scxas scy = kmpas*scyas print, name1,h1,v1,hc,vc,scx,scy,doy,tim return, a ; end ; pro cori222,name,img,cor=cor,istep=istep,pspd=pspd,htmin=htmin,htmax=htmax, $ sspeed=sspeed,fspeed=fspeed,smax=smax,angl=angl,dhts=dhts, $ az1=az1,az2=az2,ht1=ht1,ht2=ht2,ok=ok,imgg=imgg,mggg=mggg, $ xsize=xsize,ysize=ysize,iave=iave,nimage=nimage,nsmo=nsmo ; ; ; This program takes image sets from LASCO, uses the best current technique ; to difference the data from a base, and then displays this data. Using a ; mouse cursor the researcher defines an area on the image to track. This ; area is drawn on the next image and is also defined with outward motion on ; this image. The program figures out how many different preferred speeds to ; try within the limits of sspeed and fspeed and then trys them on a single ; image set before continuing. The best three answers within the given area ; are averaged together to determine the correct answer. The structures are ; NOT tracked. The program redefines the same area on the next image and ; speeds within these areas are redrawn. The LASCO data is filtered using ; the best available filter, and the following images within that same sector ; in the data set are, too. A two-D correlation program is used to follow the ; best correlation in that box outward from that area in following images. ; common DEVICES, PLOTDEV PLOTDEV = 'X' if n_elements(name) eq 0 then name = 'jan15_' ; Jan. 15 date default if n_elements(cor) eq 0 then cor = 3 ; C3 coronagraph default if n_elements(img) eq 0 then img = 2308 ; Image 2447 default if n_elements(istep) eq 0 then istep = 4 ; Default image increment if n_elements(pspd) eq 0 then pspd = 300.0 ; Pref. default speed if n_elements(htmin) eq 0 then htmin = 4.0 ; Lowest height if n_elements(htmax) eq 0 then htmax = 14.0 ; Highest height ; ; Size of the correlation window used (must be an odd number) ; if n_elements(xsize) eq 0 then xsize = 9 ; horizontal size if n_elements(ysize) eq 0 then ysize = xsize ; vertical size ; if n_elements(angl) eq 0 then angl = 20.0 ; max angle from radial if n_elements(sspeed) eq 0 then sspeed = 50.0 ; slowest speed allowed if n_elements(fspeed) eq 0 then fspeed = 1000.0 ; fastest speed allowed if n_elements(smax) eq 0 then smax = 0.1 ; lowest corr. allowed ; if n_elements(iave) eq 0 then iave = 5 ; filter tangential ; pixel average away ; from center ; if n_elements(dhts) eq 0 then dhts = 1.0 ; Height step multiplier if n_elements(az1) eq 0 then az1 = 50 ; bigger azimuth (1) if n_elements(az2) eq 0 then az2 = 40 ; lesser azimuth (2) if n_elements(ht1) eq 0 then ht1 = 4.0 ; smaller height (1) if n_elements(ht2) eq 0 then ht2 = 5.0 ; greater height (2) if n_elements(nimage) eq 0 then nimage = 33 ; number of image sets if n_elements(nsmo) eq 0 then nsmo = 0 ; smoothing of image rs = 695980.0 ; Sun radius in km ht1 = ht1*rs ; ht1 in km ht2 = ht2*rs ; ht2 in km htmin = htmin*rs ; htmin in km htmax = htmax*rs ; htmax in km ; if n_elements(ok) eq 0 then ok = 0 ; if 1, automatic run if n_elements(imgg) eq 0 then imgg = 1 ; if 1, images plotted if n_elements(mggg) eq 0 then mggg = 5 ; if 1, corrs. plotted iok = 0 if imgg ne 1 then iok = 1 ; badspv = -5.0 ; value in array if speed correlation is out of limit badang = 999.0 ; value in array if speed correlation is out of limit badamp = -5.0 ; value in array if speed correlation is out of limit ; dayins = 24.*60.*60. ; number of seconds in a day radtd = 57.29578 ; degrees in a radian itrigger = 0 itrigger2 = 0 ; dir1 = '/home/cedar/igse/bjackson/data/' dir2 = '/home/cedar/igse/bjackson/results/' fini1 = '.dif' fini2 = '.dat' openw, one, dir2+name+strcompress(img,/rem)+fini2,/get_lun print, ' RUN ON ', name+strcompress(img,/rem)+fini2 printf, one, ' RUN ON ', name+strcompress(img,/rem)+fini2 ; nam1 = name name = dir1+name ; maxini = 2000 ini = 0 istar = fltarr(2*xsize+1,2*ysize+1) istar(*,*) = -10000.0 if mggg eq 2 then begin restore, "idlsave.dat" ; use saved star data if mggg=2 print, 'Star restored from idlsave.dat' printf, one, 'Star restored from idlsave.dat' endif imgs1 = fltarr(maxini) imgs2 = fltarr(maxini) psps = fltarr(maxini) pangs = fltarr(maxini) pamps = fltarr(maxini) htms = fltarr(maxini) aznaz1 = fltarr(maxini) aznaz2 = fltarr(maxini) npspn = 1 pcos = 1.0 psin = 0.0 pang = 0.0 ht1s = ht1 ht2s = ht2 htdn = dhts*(ht2s - ht1s) az1s = az1 az2s = az2 psp = pspd ; start: ; ; Determine different parameters for this specific image. ; is = 0 ; par = param(name,img,is,fini1,namea,hsize1,vsize1,hcent1,vcent1,scx1,scy1, $ day1,time1) ; print, 'hcent1 vcent1 scx1 scy1',hcent1,vcent1,scx1,scy1 ; ; Plot the differenced image on a window ; ; Best guess of first image brightness ; ; May 24 if nam1 eq 'may24_' then begin mxbv1 = 611 mnbv1 = 499 endif ; ; Jan 15 if nam1 eq 'jan15_' then begin mxbv1 = 40 mnbv1 = 20 endif ; ; Feb 27 if nam1 eq 'feb27_' then begin mxbv1 = 40 mnbv1 = -10 hcent1 = 513.634 vcent1 = 0.0 endif ; ; Mar 02 if nam1 eq 'mar02_' then begin mxbv1 = 40 mnbv1 = -10 hcent1 = 513.634 vcent1 = 0.0 endif ; up: img1 = image(namea,hsize1,vsize1,mxbv1,mnbv1,imgg,a) ; ; Determine the area of the image to analyze (azimuth and height values) ; if imgg eq 1 then begin getc = getcp(hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,ok) if ok eq 0 then begin ht1s = ht1 ; Beginning height 1 ht2s = ht2 ; Beginning height 2 az1s = az1 ; Beginning azimuth 1 az2s = az2 ; Beginning azimuth 2 endif endif ; naz = abs(az1 - az2) naz = fix(naz) az1naz = fltarr(naz+1) az2naz = fltarr(naz+1) bestsp0 = fltarr(naz+1) bestsp1 = fltarr(naz+1) bestsp2 = fltarr(naz+1) bestamp0 = fltarr(naz+1) bestamp1 = fltarr(naz+1) bestamp2 = fltarr(naz+1) bestang0 = fltarr(naz+1) bestang1 = fltarr(naz+1) bestang2 = fltarr(naz+1) for i=0,naz do begin iaz = fix(az1 + 0.5) - i az1naz(i) = float(iaz + 0.5) az2naz(i) = float(iaz - 0.5) endfor ; ; Determine the area of the second image to analyze. ; ; Plot the first differenced image on a window with the area ; superimposed. ; dtime = 0.0 ; fd = findmx(az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,xsize,ysize,iave, $ minx,maxx,miny,maxy) ; print, 'minx maxx miny maxy', minx,maxx,miny,maxy ; mxmn = mxmnbv(a,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ minx,maxx,miny,maxy,mxbv1,mnbv1) print, 'az1 az2',az1,az2 print, 'mxbv1 mnbv1', mxbv1,mnbv1 ; imga=imagea(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp) ; print, 'az1 az2',az1,az2 if imgg eq 1 then iok = getok(hsize1,vsize1,ok) if not iok then goto, up ; ; Determine the parameters for the second image. ; par = param(name,img,istep,fini1,nameb,hsize2,vsize2,hcent2,vcent2, $ scx2,scy2,day2,time2) add = 0.0 if day2 lt day1 then add = (day1-day2)*dayins if day2 gt day1 then add = (day2-day1)*dayins time2 = time2 + add ; time of image 2 in seconds dtime = time2-time1 ; time difference between images ; ; ; determine the updated height and azimuth differences to be used on ; the second image for the next correlations ; htd = psp*dtime*pcos ; ; print, htd,psp,dtime,pcos ; ; Determine the next image area to filter and analyze ; filta = abs(atan((1.5*xsize),ht1/scx1))*radtd ; extra azimuth in deg filtb1 = 1.5*xsize*scx1 ; extra -height value in km filtb2 = 1.5*xsize*scx1 ; extra +height value in km hts = sspeed*dtime*pcos htf = fspeed*dtime*pcos htm = (ht1 + ht2)/2.0 ; Mean height of first analysis azd = abs(atan((psin*hts),(htm + hts))*radtd) ; up2: ; next = nextim(az1,az2,ht1,ht2,filta,filtb1,filtb2, $ azd,hts,htf,az11,az22,ht11,ht22) print, 'az11 az22 ht11 ht22',az11,az22,ht11/rs,ht22/rs ; ; Find the min and max x and y values for the area analyzed on the second image ; ; Best guess of second image brightness ; May 24 if nam1 eq 'may24_' then begin mxbv2 = 611 mnbv2 = 499 endif ; ; Jan 15 if nam1 eq 'jan15_' then begin mxbv2 = 40 mnbv2 = 20 endif ; ; Feb 27 if nam1 eq 'feb27_' then begin mxbv2 = 40 mnbv2 = -10 hcent2 = 513.634 vcent2 = 0.0 endif ; ; Mar 02 if nam1 eq 'mar02_' then begin mxbv2 = 40 mnbv2 = -10 hcent2 = 513.634 vcent2 = 0.0 endif ; fd=findmx(az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2,xsize,ysize,iave, $ minx11,maxx11,miny11,maxy11) print, minx11,maxx11,miny11,maxy11 ; spcw = xsize*scx1/dtime ; ; Plot the second differenced image on a window with the new area ; superimposed. ; img2 = image(nameb,hsize2,vsize2,mxbv2,mnbv2,imgg,b) ; mxmn = mxmnbv(b,hsize2,vsize2,az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2, $ minx11,maxx11,miny11,maxy11,mxbv2,mnbv2) ; ; locate the updated height and the azimuth location on the second image ; ht1d = ht1 + htd ht2d = ht2 + htd az1d = az1 + azd az2d = az2 + azd ; imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az1d,az2d,ht1d,ht2d,imgg, $ xsize,ysize,az1,ht2,time1,time2,dtime,psp) imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time2,dtime,psp) imga=imagea(hsize2,vsize2,hcent2,vcent2,scx2,scy2,az11,az22,ht11,ht22,imgg, $ xsize,ysize,az1,ht2,time1,time2,dtime,psp) print, 'az1 az2',az1,az2 print, 'az11 az22',az11,az22 ; print, 'Speed and angle on the image',psp,pang ; if imgg eq 1 then iok = getok(hsize2,vsize2,ok) if not iok and imgg eq 1 then begin gsp = getspeed(hcent2,vcent2,scx2,scy2,ht1,ht2,dtime,psp,ok) goto, up2 endif ok = 1 ; ; Filter the first and second image areas. The area beyond this region is ; not filtered. ; ; First image ; smo = smoo(a,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,ar) if nsmo ge 3 then begin ar = smooth(ar,nsmo) print, 'Image 1 smoothed',nsmo,' by',nsmo printf, one, 'Image 1 smoothed',nsmo,' by',nsmo endif ft1 = filtr(ar,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,i1) ; ; Second image ; smo = smoo(b,hsize2,vsize2,xsize,ysize,hcent2,vcent2,scx2,scy2, $ az11,az22,ht11,ht22,minx11,maxx11,miny11,maxy11,br) if nsmo ge 3 then begin br = smooth(br,nsmo) print, 'Image 2 smoothed',nsmo,' by',nsmo printf, one, 'Image 2 smoothed',nsmo,' by',nsmo endif ; ft2 = filtr(br,hsize2,vsize2,iave,xsize,ysize,hcent2,vcent2,scx2,scy2, $ az11,az22,ht11,ht22,minx11,maxx11,miny11,maxy11,i2) ; ; Correlate from one location to the other in the two images using the ; nominal speed value. ; nkx = long(maxx-minx+1) nky = long(maxy-miny+1) nk = nkx*nky print, 'nk', nk spv = fltarr(hsize1,vsize1) ampv = fltarr(hsize1,vsize1) angv = fltarr(hsize1,vsize1) spv(*,*) = sspeed ampv(*,*) = smax angv(*,*) = 0.0 pspeedxy = fltarr(nk) xcent = fltarr(nk) ycent = fltarr(nk) i11s = fltarr(2*xsize+1,2*ysize+1) i22s = fltarr(2*xsize+1,2*ysize+1) ccfs = fltarr(xsize,ysize) amaxs = -1.0 ; pspeedxy(*) = badspv ; hh1 = time1/(60.*60.) hh1 = fix(hh1) mm1f = time1/60. - hh1*60. mm1 = fix(mm1f) ss1 = mm1f*60. - mm1*60. hh2 = time2/(60.*60.) hh2 = fix(hh2) mm2f = time2/60. - hh2*60. mm2 = fix(mm2f) ss2 = mm2f*60. - mm2*60. print, ' ' print, 'first image is at ', hh1, mm1, ss1 printf, one, 'first image is at', hh1, mm1, ss1 print, 'second image is at', hh2, mm2, ss2 printf, one, 'second image is at', hh2, mm2, ss2 print, 'correlation window width is equivalent to ',spcw,' km/s' print, 'azimuths ', az1, az2,' hts ', ht1/rs, ht2/rs printf, one, 'correlation window width is equivalent to ',spcw,' km/s' printf, one, 'azimuths', az1, az2,' hts', ht1/rs, ht2/rs ; bestsp0(*) = 0.0 bestsp1(*) = 0.0 bestsp2(*) = 0.0 bestamp0(*) = 0.0 bestamp1(*) = 0.0 bestamp2(*) = 0.0 bestang0(*) = 0.0 bestang1(*) = 0.0 bestang2(*) = 0.0 ; ; Find the best preferred speed values to use to propagate to the next ; image (psp is cycled from sspeed to fspeed at distances of 1/3 xsize). ; sxsize3 = spcw/3.0 ; 1/3 the speed equivalent to the length of xsize psp = sspeed again: print, ' psp', psp ; cx = cimage(i1,i2,hsize1,vsize1,hsize2,vsize2,hcent1,vcent1,hcent2,vcent2, $ scx1,scy1,scx2,scy2,psp,dtime,az1,az2,ht1,ht2, $ minx,maxx,miny,maxy,amaxs,i11s,i22s,ccfs,istar, $ pspeedxy,xcent,ycent,kk,itrigger,imgg,mggg, $ az11,az22,ht11,ht22,xsize,ysize,one, $ angl,sspeed,fspeed,smax,badspv,badang,badamp, $ spv,angv,ampv) ; ; Determine the highest amplitude correlations of speed and angle values ; within the area defined by az1naz, az2naz and ht1, ht2. ; for i=0,naz do begin bests0 = bestsp0(i) bests1 = bestsp1(i) bests2 = bestsp2(i) bestam0 = bestamp0(i) bestam1 = bestamp1(i) bestam2 = bestamp2(i) bestan0 = bestang0(i) bestan1 = bestang1(i) bestan2 = bestang2(i) ; fd = findmx(az1naz(i),az2naz(i),ht1,ht2,hcent1,vcent1,scx1,scy1, $ xsize,ysize,iave,mnx,mxx,mny,mxy) ; bv=bestv(spv,ampv,angv,hsize1,vsize1,hcent1,vcent1,scx1,scy1,xsize,ysize, $ az1naz(i),az2naz(i),ht1,ht2,mnx,mxx,mny,mxy, $ bests0,bests1,bests2,bestam0,bestam1,bestam2,bestan0,bestan1,bestan2) ; ; print, $ ;az1naz(i),az2naz(i),ht1/rs,ht2/rs,bestam0,bests0,bestan0,bestam1,bests1,bestan1 bestsp0(i) = bests0 bestsp1(i) = bests1 bestsp2(i) = bests2 bestamp0(i) = bestam0 bestamp1(i) = bestam1 bestamp2(i) = bestam2 bestang0(i) = bestan0 bestang1(i) = bestan1 bestang2(i) = bestan2 endfor ; if psp gt fspeed then goto, next psp = psp + sxsize3 goto, again next: psp = pspd ; print, 'az1naz az2naz', az1naz, az2naz printf, one, 'az1naz az2naz', az1naz, az2naz print, 'speed0 =', bestsp0 printf, one, 'speed0 =', bestsp0 print, 'speed1 =', bestsp1 printf, one, 'speed1 =', bestsp1 print, 'speed2 =', bestsp2 printf, one, 'speed2 =', bestsp2 ; print, 'amplitude0 =', bestamp0 printf, one, 'amplitude0 =', bestamp0 print, 'amplitude1 =', bestamp1 printf, one, 'amplitude1 =', bestamp1 print, 'amplitude2 =', bestamp2 printf, one, 'amplitude2 =', bestamp2 ; print, 'angle0 =', bestang0 printf, one, 'angle0 =', bestang0 print, 'angle1 =', bestang1 printf, one, 'angle1 =', bestang1 print, 'angle2 =', bestang2 printf, one, 'angle2 =', bestang2 ; ; Plot the best speed values within the correlated area superposed on the ; image ; if imgg eq 1 or imgg eq 2 then begin if mggg ge 1 then begin MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; mincs = Blueish(0) maxcs = Blueish(1) mincg = Reddish(0) maxcg = Greenish(1) mincp = Reddish(0) maxcp = Reddish(1) endif endif ; img10 = image10(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ imgg,minx,maxx,miny,maxy,spv,sspeed,fspeed,mxbv1,mnbv1,mincs,maxcs) ; ; Overplot the boundary and the speed color scale box ; imgb = imageb(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp, $ sspeed,fspeed,mincs,maxcs) ; ; Overplot the speed scale on the box ; imgc = imagec(hsize1,vsize1,sspeed,fspeed,mincs,maxcs,imgg) ; ; Plot the best amplitude values within the correlated area superposed on the ; image ; img10 = image10(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,imgg, $ minx,maxx,miny,maxy,ampv,0.0,1.0,mxbv1,mnbv1,mincp,maxcp) ; ; Overplot the boundary and the amplitude color scale box ; imgb = imageb(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp, $ 0.0,1.0,mincp,maxcp) ; ; Overplot the amplitude scale on the box ; imgc = imagec(hsize1,vsize1,0.0,1.0,mincp,maxcp,imgg) ; ; Plot the best angle values within the correlated area superposed on the ; image ; img10 = image10(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1,imgg, $ minx,maxx,miny,maxy,angv,-angl,angl,mxbv1,mnbv1,mincg,maxcg) ; ; Overplot the boundary and the angle color scale box ; imgb = imageb(hsize1,vsize1,hcent1,vcent1,scx1,scy1,az1,az2,ht1,ht2,imgg, $ xsize,ysize,az1,ht2,time1,time1,dtime,psp, $ -angl,angl,mincg,maxcg) ; ; Overplot the angle scale on the box ; imgc = imagec(hsize1,vsize1,-angl,angl,mincg,maxcg,imgg) ; print, 'npspn nimage', npspn, nimage ; ; Plot out the best correlation of the bunch ; ; *** ; if imgg le 2 then begin if mggg ge 1 and itrigger le 2 then begin if ht1 gt 7.5 then begin print, 'set of images correlated and correlation', amaxs printf, one, 'set of images correlated and correlation', amaxs ; print, 'Contouring first image' WINDOW, /free, xsize = 420, ysize=420 mini1 = min(i11s(*,*)) maxi1 = max(i11s(*,*)) print, '1st image mini1 maxi1', mini1, maxi1 printf, one, '1st image mini1 maxi1', mini1, maxi1 i11s(*,*) = i11s(*,*) - mini1 i11s(*,*) = i11s(*,*)*(0.5)/(maxi1-mini1) contour, i11s, levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [81,84,88,92,96,100,104,108,112,116,120], /fill, $ xstyle=1, ystyle=1, pos=[0.1,0.1,0.9,0.9] ; print, 'Contouring second image' WINDOW, /free, xsize = 420, ysize=420 mini2 = min(i22s(*,*)) maxi2 = max(i22s(*,*)) print,'2nd image mini2 maxi2', mini2, maxi2 printf, one, '2nd image mini2 maxi2', mini2, maxi2 i22s(*,*) = i22s(*,*) - mini2 i22s(*,*) = i22s(*,*)*(0.5)/(maxi2-mini2) contour, i22s, levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [81,84,88,92,96,100,104,108,112,116,120], /fill, $ xstyle=1, ystyle=1, pos=[0.1,0.1,0.9,0.9] ; print, 'Contouring correlation height for a sample' WINDOW, /free, xsize = 186, ysize=186 minic = min(ccfs(*,*)) maxic = max(ccfs(*,*)) print,'Correlation minic maxic', minic, maxic printf, one, 'Correlation minic maxic', minic, maxic cc = fltarr(xsize,ysize) cc(*,*) = ccfs(*,*) - minic cc(*,*) = cc(*,*)*(0.5)/(maxic-minic) contour, cc, levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [121,124,128,132,136,140,144,148,152,156,160], /fill, $ xstyle=1, ystyle=1, pos=[0.1,0.1,0.9,0.9] ; itrigger = itrigger + 1 endif endif endif ; ; *** ; ; Save a stellar image if there is one. ; if mggg eq 3 then begin istar(*,*) = i22s(*,*) ; copy every second image to file save, istar print, 'Every second image copied to idlsave.dat file' printf, one, 'Every second image copied to idlsave.dat file' endif if itrigger2 eq 0 and mggg eq 1 then begin istar(*,*) = i22s(*,*) ; copy the stellar image to the istar array save, istar ; save only the first stellar image to file print, 'Only the first second image is copied to idlsave.dat file' printf, one, 'Only the first second image is copied to idlsave.dat file' endif ; itrigger2 = itrigger2 + 1 ; ; Do the correlations once more if npspn lt nimage ; ; for i=0,naz do begin pspi = (bestsp0(i)*bestamp0(i)+bestsp1(i)*bestamp1(i)+bestsp2(i)*bestamp2(i))/ $ (bestamp0(i) + bestamp1(i) + bestamp2(i)) pang=(bestang0(i)*bestamp0(i)+bestang1(i)*bestamp1(i)+bestang2(i)*bestamp2(i))/$ (bestamp0(i) + bestamp1(i) + bestamp2(i)) pamp=(bestamp0(i)*bestamp0(i)+bestamp1(i)*bestamp1(i)+bestamp2(i)*bestamp2(i))/$ (bestamp0(i) + bestamp1(i) + bestamp2(i)) ; ; ; determine the updated height and azimuths to be used on the first image ; of the next set for the next correlations ; ; pcosi = cos(pang/radtd) psini = sin(pang/radtd) if pcosi ne 0 then begin htd = pspi*dtime*pcosi endif ; htmms = htm + htd ; Mean height from speeds on images ; ; azd = atan((psini*htdn),(htm + htdn))*radtd ; az1 = az1 + azd ; New starting azimuth 1 az2 = az2 + azd ; New starting azimuth 2 ; print, 'I1 best speed',img,' at',pamp,' at',pspi,' at angle', pang,' at', $ ; htmms/rs,' Rs' ; printf, one, 'I1 best speed',img,' at',pamp,' at',pspi,' at angle', pang, $ ; ' at', htmms/rs,' Rs' imgs1(ini) = img imgs2(ini) = img + istep psps(ini) = pspi pamps(ini) = pamp pangs(ini) = pang htms(ini) = htmms/rs imgs1(ini) = img aznaz1(ini) = az1naz(i) aznaz2(ini) = az2naz(i) if ini ge maxini then goto, end1 ini = ini + 1 endfor ; ht1 = ht1 + htdn ; New starting height 1 ht2 = ht2 + htdn ; New starting height 2 ; if ht1 gt htmax then begin for i=0,ini-1 do begin ; just in case, print summary print, format='(f5.0,2f6.1,f7.2,f7.3,f7.4,f7.3)', $ imgs1(i),aznaz1(i),aznaz2(i),psps(i),pangs(i),pamps(i),htms(i) print, one, format='(f5.0,2f6.1,f7.2,f7.3,f7.4,f7.3)', $ imgs1(i),aznaz1(i),aznaz2(i),psps(i),pangs(i),pamps(i),htms(i) endfor npspn = npspn + 1 ; Bump image counter img = img + istep ; New starting image ht1 = ht1s ; Beginning height 1 ht2 = ht2s ; Beginning height 2 az1 = az1s ; Beginning azimuth 1 az2 = az2s ; Beginning azimuth 2 endif ; if npspn gt nimage then goto, end1 ; last time through print, 'New image1 az1 az2 ht1 ht2 psp pang',img,az1,az2,ht1/rs,ht2/rs,psp,pang goto, start ; end1: print, 'Image and speed summary' printf, one, 'Image and speed summary' for i=0,ini-1 do begin print, format='(f5.0,2f6.1,f7.2,f7.3,f7.4,f7.3)', $ imgs1(i),aznaz1(i),aznaz2(i),psps(i),pangs(i),pamps(i),htms(i) print, one, format='(f5.0,2f6.1,f7.2,f7.3,f7.4,f7.3)', $ imgs1(i),aznaz1(i),aznaz2(i),psps(i),pangs(i),pamps(i),htms(i) endfor ; free_lun, one ; end