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 bestvn,sp,amp,ang,nx,ny,xc,yc,scx,scy,xsize,ysize, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestsp0,bestamp0,bestang0 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 ;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) endif endif endif ENDFOR ENDFOR ;print, az1,az2,ht1,ht2,minx,maxx,miny,maxy,bestamp0,bestsp0,bestang0 ; end ; function ccfss, ccfs,xxi,yyi,iazv ; ; This function interpolates over the pixels near the speed location on the ; correlation function. ; xxif = fix(xxi) yyif = fix(yyi) xxd = xxi - xxif yyd = yyi - yyif ; cca = (ccfs(xxif+1,yyif) - ccfs(xxif,yyif))*xxd + $ ccfs(xxif,yyif) ccb = (ccfs(xxif+1,yyif+1) - ccfs(xxif,yyif+1))*xxd + $ ccfs(xxif,yyif+1) ccfss = (ccb - cca)*yyd + cca ; return, ccfss ; end ; function sumv, spv,ampv,angv,i1,i2,mm,mm0,ccfst,nccfst,spcw,spbin,xxis,yyjs, $ nx,ny,xc,yc,iazv,scx,scy,xsize,ysize, $ az1,az2,iaz,ht1,ht2,mnx,mxx,mny,mxy,spsum,angsum, $ dif,difv,dif0,ndifv,nsum,nspmax,nspmin,iflg,psp, $ angl,sspeed,fspeed,smax ; radtd = 57.29578 sperpix = spcw/(xsize-1.0) azm = (az1 + az2)/2.0 azmr = azm/radtd dyyj = -spbin*sin(azmr)/sperpix dxxi = -spbin*cos(azmr)/sperpix nmid = (xsize-1.0)/2.0 ximlim = xsize/6 yimlim = ysize/6 xiplim = xsize - xsize/6 - 1.0 yiplim = ysize - ysize/6 - 1.0 xxi = 0.0 yyj = 0.0 dif = 0.0 dif0 = 0.0 difv = 0.0 ndifv = 0.0 xsized6 = xsize/6 ispdmx = fix(psp/spbin + 0.5) ; ccfstn = fltarr(xsize,ysize) ; nh1m = nx - 1 nv1m = ny - 1 ; ivv = ysize if mny gt ivv then ivv = mny ; ivp = nv1m - ysize if mxy lt ivp then ivp = mxy ; ihh = xsize if mnx gt ihh then ihh = mnx ; ihp = nh1m - xsize if mxx lt ihp then ihp = mxx ; ; ivv = mny ; ivp = mxy ; ihh = mnx ; ihp = mxx ; 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 ;print, az1,az2,ht1,ht2,mnx,mxx,mny,mxy,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 ; ; iflg = 1 is the nominal where the correlation maximum is used to ; determine the speed. Otherwise the center of the correlation is used. ; if iflg eq 1 then xxi = xxis(ixx,iyy) if iflg eq 1 then yyj = yyjs(ixx,iyy) ; if xxis(ixx,iyy) ne -999.0 and $ abs(angv(ixx,iyy)) le angl then begin ; ; Determine image difference sums of the successfully correlated images ; difv = difv + (i1(ixx,iyy) + i2(ixx,iyy))/2.0 dif0 = dif0 + mm0(ixx,iyy) dif = dif + mm(ixx,iyy) ndifv = ndifv + 1.0 ; ; print, xxis(ixx,iyy),yyjs(ixx,iyy),angv(ixx,iyy),spv(ixx,iyy),ampv(ixx,iyy) ; Got to here ; nccf = fix(nccfst(ixx,iyy)) ccfstn(*,*) = ccfst(*,*,nccf) dang = angv(ixx,iyy) if iflg eq 1 then dyyj = -spbin*sin((dang+azm)/radtd)/sperpix if iflg eq 1 then dxxi = -spbin*cos((dang+azm)/radtd)/sperpix spdmax = spv(ixx,iyy) xxim = xxi + nmid yyim = yyj + nmid if iflg eq 1 then ispdmx = fix(spdmax/spbin + 0.5) if ispdmx ge 0 and ispdmx le nspmax then begin if xxim gt ximlim and xxim lt xiplim then begin if yyim gt yimlim and yyim lt yiplim then begin spsum(ispdmx) = spsum(ispdmx) + $ ccfss(ccfstn,xxim,yyim) ; ampv(ixx,iyy) angsum(ispdmx) = angsum(ispdmx) + angv(ixx,iyy) nsum(ispdmx) = nsum(ispdmx) + 1.0 endif endif endif ; ; xxx=0 ; if xxx ne 0 then begin FOR i=1,xsized6 do BEGIN xipm = xxi + nmid + dxxi*i ximm = xxi + nmid - dxxi*i yjpm = yyj + nmid + dyyj*i yjmm = yyj + nmid - dyyj*i ip1 = ispdmx + i im1 = ispdmx - i if ip1 le nspmax then begin if xipm gt ximlim and xipm lt xiplim then begin if yjpm gt yimlim and yjpm lt yiplim then begin spsum(ip1) = spsum(ip1) + $ ccfss(ccfstn,xipm,yjpm) angsum(ip1) = angsum(ip1) + angv(ixx,iyy) nsum(ip1) = nsum(ip1) + 1.0 endif endif endif if im1 ge nspmin and im1 le nspmax then begin if ximm gt ximlim and ximm lt xiplim then begin if yjmm gt yimlim and yjmm lt yiplim then begin spsum(im1) = spsum(im1) + $ ccfss(ccfstn,ximm,yjmm) angsum(im1) = angsum(im1) + angv(ixx,iyy) nsum(im1) = nsum(im1) + 1.0 endif endif endif ENDFOR ; endif endif endif endif ENDFOR ENDFOR ; 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,ilnh,ilnv, $ ieh,iev,ccfst,nccfst,amaxb,mggg,badsum, $ az11,az22,ht11,ht22,htbb1,htbb2,xxiss,yyjss,xsize,one, $ angl,sspeed,fspeed,smax,badspeed,badangle,badcorr, $ cspeed,cangl,cmax,xxis,yyjs,msts,mstl,mste ; ; This function obtains sets of image arrays of the proper size and ; displacement to be used in a correlation, and then does the correlation. ; If a filtered stellar image has been provided (istar(0,0,naz) ne -10000.0) ; then the stellar image is compared with both the first and second images ; correlated. If a horizontal line has been provided ; (ilnh(0,0,naz) ne -10000.0) then the line is compared with both the first ; and second images correlated. If a vertical line has been provided ; (ilnh(0,0,naz) ne -10000.0) then the line is compared with both the first ; and second images correlated. ; radtd = 57.29578 ; rs = 695980.0 ; Sun radius in km goodcorr = 1.0 xys = badspeed dang = badangle amax = badcorr ; ysize = xsize 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 ; nkx = long(maxx-minx+1) nky = long(maxy-miny+1) nk = nkx*nky ; 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 st1dmaxd = 100000. st2dmaxd = st1dmaxd print, 'i1max, i1min, st1dmaxd', i1max, i1min, st1dmaxd ; nxstar = 0L nxlnh = 0L nxedgh = 0L nxbounds = 0L nxsame = 0l nxspeed = 0L nxang = 0L nxamp = 0L nxxx = 0L nxxxx = 0L ncount = 0L kk = 0L FOR iy=ivvv,ivvp,jstep do BEGIN FOR ix=ihhh,ihhp,istep do BEGIN if amaxb(ix,iy) eq goodcorr then begin ; ; Determine the boundary within which to do the correlation on image 1 ; bd1 = bound(ix,iy,xc1,yc1,scx1,scy1,azbb,htbb) iazmth = fix(az1 - azbb + 0.5) 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. ; ncount = ncount + 1 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 ; iflag0 = 1 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). ; iflag7 = 0 iflag6 = 0 iflag5 = 0 if iflag0 eq 1 then begin if iflag1 eq 0 or iflag2 eq 0 or iflag3 eq 0 then nxsame = nxsame + 1 if iflag1 eq 1 then begin if iflag2 eq 1 then begin if iflag3 eq 1 then begin if iflag4 eq 1 then begin ; ; Assume values good unless otherwise flagged or known ; xys = sspeed + 0.01 dang = angl - angl/2.0 amax = smax + 0.01 ; 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 iflag0 = 0 nxbounds = nxbounds + 1 endif else begin ; ; Now check to see if the images are starlike ; ; print, 'istar(0,0,iazmth)', istar(0,0,iazmth) ; ; msts = 0.55 ; If corr. higher than this, the image contains a star amax1s = 0.0 amax2s = 0.0 if istar(0,0,iazmth) 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(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax1s) endif else begin amax1s = 1.0 endelse ; if abs(amax1s) gt msts then print, 'Image 1 contains a star - amax1s msts', amax1s, msts ; ; 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(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax2s) endif else begin amax2s = 1.0 endelse ; if abs(amax2s) gt msts then print, 'Image 2 contains a star - amax2s msts', amax2s, msts endif ; if abs(amax1s) le msts and abs(amax2s) le msts then iflag5 = 1 ; image not a star ; if iflag5 eq 0 then begin nxstar = nxstar + 1 goto, endcheck endif ; ; Now check to see if the image contains a horizontal or vertical line ; ; mstl = 0.25 ; If corr. higher than this, the image contains a line amax1s = 0.0 amax2s = 0.0 if ilnh(0,0,iazmth) ne -10000.0 then begin ; ; Compare image 1 with a horizontal line ; i1Xst = xyoff(i11,ilnh(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax1s) ; if abs(amax1s) gt mstl then print, 'Image 1 contains a horizontal line - amax1s mstl', amax1s, mstl ; ; Compare image 2 with a horizontal line ; i2Xst = xyoff(i22,ilnh(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax2s) ; if abs(amax2s) gt mstl then print, 'Image 2 contains a horizontal line - amax2s mstl', amax2s, mstl endif ; if abs(amax1s) le mstl and abs(amax2s) le mstl then iflag6 = 1 ; image not a line ; if iflag6 eq 0 then begin nxlnh = nxlnh + 1 goto, endcheck endif ; amax1s = 0.0 amax2s = 0.0 if ilnv(0,0,iazmth) ne -10000.0 then begin ; ; Compare image 1 with a vertical line ; i1Xst = xyoff(i11,ilnv(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax1s) ; if abs(amax1s) gt mstl then print, 'Image 1 contains a vertical line - amax1s mstl', amax1s, mstl ; ; Compare image 2 with a vertical line ; i2Xst = xyoff(i22,ilnv(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax2s) ; if abs(amax2s) gt mstl then print, 'Image 2 contains a vertical line - amax2s mstl', amax2s, mstl endif ; if abs(amax1s) le mstl and abs(amax2s) le mstl then iflag6 = 1 ; image not a line ; if iflag6 eq 0 then begin nxlnh = nxlnh + 1 goto, endcheck endif ; ; ; There is no star or line contained in the image, check for an edge ; ; mste = 0.25 ; If corr. higher than this, the image contains a line amax1s = 0.0 amax2s = 0.0 if ieh(0,0,iazmth) ne -10000.0 then begin ; ; Compare image 1 with a horizontal edge ; i1Xst = xyoff(i11,ieh(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax1s) ; if amax1s gt mste then print, 'Image 1 contains a horizontal edge - amax1s mste', amax1s, mste ; ; Compare image 2 with a horizontal edge ; i2Xst = xyoff(i22,ieh(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax2s) ; if amax2s gt mste then print, 'Image 2 contains a horizontal edge - amax2s mste', amax2s, mste endif ; if abs(amax1s) le mste and abs(amax2s) le mste then iflag7 = 1 ; image not an edge ; if iflag7 eq 0 then begin nxedgh = nxedgh + 1 goto, endcheck endif ; amax1s = 0.0 amax2s = 0.0 if iev(0,0,iazmth) ne -10000.0 then begin ; ; Compare image 1 with a vertical edge ; i1Xst = xyoff(i11,iev(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax1s) ; if abs(amax1s) gt mste then print, 'Image 1 contains a vertical edge - amax1s mste', amax1s, mste ; ; Compare image 2 with a vertical edge ; i2Xst = xyoff(i22,iev(*,*,iazmth),xsize,ysize,xsize,ysize,xxii, $ yyjj,amax2s) ; if abs(amax2s) gt mste then print, 'Image 2 contains a vertical edge - amax2s mste', amax2s, mste endif ; if abs(amax1s) le mste and abs(amax2s) le mste then iflag7 = 1 ; image not a line ; if iflag7 eq 0 then begin nxedgh = nxedgh + 1 goto, endcheck endif ; endcheck: ; 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 ; xys9 = xys dang9 = dang amax9 = amax if iflag0 eq 0 then begin xys = badspeed dang = badangle amax = badcorr 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 abs(xys) lt 100.0 and amax gt 0.5 and pspeed eq 0.0 then begin amaxb(ix,iy) = 0.0 ncount = ncount - 1 endif ; if iflag5 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag6 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif if iflag7 eq 0 then begin xys = badspeed dang = badangle amax = badcorr endif ; if xys ne badspeed then begin if abs(dang) gt angl then begin nxang = nxang + 1 xys = badspeed dang = badangle amax = badcorr endif if amax lt smax then begin nxamp = nxamp + 1 xys = badspeed dang = badangle amax = badcorr endif if xys lt sspeed then begin nxspeed = nxspeed + 1 xys = badspeed dang = badangle amax = badcorr endif if xys gt fspeed then begin nxspeed = nxspeed + 1 xys = badspeed dang = badangle amax = badcorr endif endif endif endif endif endif endif ; ; Save the maximum amplitude images arrays and the correlation array plus the ; non-bad values of all the speed, angle and amplitude for each image pixel ; center correlated. Also save the cordinates of the pixel centers. ; ; if iflag5 eq 1 and iflag6 eq 1 and iflag7 eq 1 then $ ; print, iazmth,dang9,amax9,xys9,iflag0,iflag1,iflag2,iflag3,iflag4,iflag5,iflag6,iflag7 ; if iflag5 eq 1 and iflag6 eq 1 and iflag7 eq 1 then $ ; printf, one, iazmth,dang9,amax9,xys9,iflag0,iflag1,iflag2,iflag3,iflag4,iflag5,iflag6,iflag7 ; if amaxb(ix,iy) ne 0.0 then begin if pspeed eq 0.0 then badsum(iazmth) = 1.0 if xys ne badspeed and dang ne badangle and amax ne badcorr then begin xxis(ix,iy) = xxi yyjs(ix,iy) = yyj ; nxxxx = nxxxx + 1L ; nccfst(ix,iy) = kk + 0.5 ccfst(*,*,kk) = ccf(*,*) ; if kk ge nk then begin print, ' kk limit exceeded' printf, one, ' kk limit exceeded' goto, endkk endif ; kk = kk + 1L ; if amax gt amaxs(iazmth) then begin amaxs(iazmth) = amax i11s(*,*,iazmth) = i11(*,*) i22s(*,*,iazmth) = i22(*,*) ccfs(*,*,iazmth) = ccf(*,*) htbb1(iazmth) = htbb htbb2(iazmth) = htb xxiss(iazmth) = xxi yyjss(iazmth) = yyj ; print, 'Max amp image- iazmth amax1s amax2s amaxs(iazmth) xys dang azbb htb htbb xxi yyj', $ iazmth, amax1s, amax2s, amaxs(iazmth), xys, dang, azbb, htb/rs, htbb/rs, xxi, yyj endif ; endif endif cspeed(ix,iy) = xys cangl(ix,iy) = dang cmax(ix,iy) = amax ; ; endif endif endif ENDFOR ENDFOR ; endkk: ; 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' print, ' ', nxlnh, ' pixels were line correlations' printf, one, ' ', nxlnh, ' pixels were line correlations' print, ' ', nxedgh, ' pixels were edge correlations' printf, one, ' ', nxedgh, ' pixels were edge correlations' print, ' ', nxsame, ' images were the same or bad' printf, one, ' ', nxsame, ' images were the same or bad' print, ' ', nxbounds, ' pixels correlated out of bounds' printf, one, ' ', nxbounds, ' pixels correlated out of bounds' print, ' ', nxspeed, ' pixels failed speed criteria' printf, one, ' ', nxspeed, ' pixels failed speed criteria' print, ' ', nxang, ' pixels failed angle criteria' printf, one, ' ', nxang, ' pixels failed angle criteria' print, ' ', nxamp, ' pixels failed amplitude criteria' printf, one, ' ', nxamp, ' pixels failed amplitude criteria' print, ' ', ncount, ' pixels passed 0 speed criteria' printf, one, ' ', ncount, ' pixels passed 0 speed criteria' ; 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 starput, stari,hsize1,vsize1,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2 ; ; This function places stellar images in the image stari in order that ; they cam be filtered and smoothed just like the real differenced ; coronagraph images ; radtd = 57.29578 ; ; Determine the position of each stellar image ; azi = 9.0 azp = azi/2.0 htm = (ht1 + ht2)/2.0 FOR az=az2+azp,az1,azi do BEGIN si = sin((az)/radtd) co = cos((az)/radtd) ix = fix(hcent1 + (htm/scx1)*co) iy = fix(vcent1 + (htm/scy1)*si) stari(ix,iy) = 2.0 stari(ix+1,iy) = 2.0 stari(ix,iy+1) = 2.0 stari(ix+1,iy+1) = 2.0 ENDFOR ; END ; function hlineput, hlnri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2 ; ; This function places horizontal line images in the image stari in order ; that they cam be filtered and smoothed just like the real differenced ; coronagraph images ; radtd = 57.29578 ; ; Determine the position of each line image ; azi = 9.0 azp = azi/2.0 htm = (ht1 + ht2)/2.0 FOR az=az2+azp,az1,azi do BEGIN si = sin((az)/radtd) co = cos((az)/radtd) ix = fix(hcent1 + (htm/scx1)*co) iy = fix(vcent1 + (htm/scy1)*si) FOR i=0, xsize+1 do BEGIN hlnri(ix+i,iy) = 2.0 hlnri(ix-i,iy) = 2.0 ENDFOR ENDFOR ; END ; function vlineput, vlnri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2 ; ; This function places vertical line images in the image stari in order ; that they cam be filtered and smoothed just like the real differenced ; coronagraph images ; radtd = 57.29578 ; ; Determine the position of each line image ; azi = 9.0 azp = azi/2.0 htm = (ht1 + ht2)/2.0 FOR az=az2+azp,az1,azi do BEGIN si = sin((az)/radtd) co = cos((az)/radtd) ix = fix(hcent1 + (htm/scx1)*co) iy = fix(vcent1 + (htm/scy1)*si) FOR j=0, ysize+1 do BEGIN vlnri(ix,iy+j) = 2.0 vlnri(ix,iy-j) = 2.0 ENDFOR ENDFOR ; END ; function hedgeput, hedgri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2 ; ; This function places horizontal edge images in the image hedgri in order ; that they cam be filtered and smoothed just like the real differenced ; coronagraph images ; radtd = 57.29578 ; ; Determine the position of each image edge ; azi = 9.0 azp = azi/2.0 htm = (ht1 + ht2)/2.0 FOR az=az2+azp,az1,azi do BEGIN si = sin((az)/radtd) co = cos((az)/radtd) ix = fix(hcent1 + (htm/scx1)*co) iy = fix(vcent1 + (htm/scy1)*si) FOR j=0,ysize do BEGIN FOR i=0, xsize do BEGIN hedgri(ix+i,iy+j) = 2.0 hedgri(ix-i,iy+j) = 2.0 ENDFOR ENDFOR ENDFOR ; END ; function vedgeput, vedgri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2 ; ; This function places vertical edge images in the image vedgri in order ; that they cam be filtered and smoothed just like the real differenced ; coronagraph images ; radtd = 57.29578 ; ; Determine the position of each image edge ; azi = 9.0 azp = azi/2.0 htm = (ht1 + ht2)/2.0 FOR az=az2+azp,az1,azi do BEGIN si = sin((az)/radtd) co = cos((az)/radtd) ix = fix(hcent1 + (htm/scx1)*co) iy = fix(vcent1 + (htm/scy1)*si) FOR i=0,xsize do BEGIN FOR j=0, ysize do BEGIN vedgri(ix+i,iy+j) = 2.0 vedgri(ix+i,iy-j) = 2.0 ENDFOR ENDFOR ENDFOR ; END ; function extrims, starims,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,istar ; ; This function extracts stellar images from the filtered and smoothed ; image starims ; radtd = 57.29578 ; azi = 9.0 iazd2 = fix(azi/2.0) iazi = fix(azi + 0.1) azp = azi/2.0 naz = fix(az1 - azp - az2) htm = (ht1 + ht2)/2.0 FOR az=az2+azp,az1,azi do BEGIN si = sin((az)/radtd) co = cos((az)/radtd) ; ; Center of the stellar image ; ix = fix(hcent1 + (htm/scx1)*co) iy = fix(vcent1 + (htm/scy1)*si) ; ivv = iy-ysize ivp = iy+ysize ihh = ix-xsize ihp = ix+xsize ; iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN if naz ge 0 then istar(ixh,iyv,naz) = starims(ixx,iyy) if naz+iazd2+1 ge 0 then istar(ixh,iyv,naz+iazd2+1) = starims(ixx,iyy) FOR i=1,iazd2 do BEGIN if naz+i ge 0 then istar(ixh,iyv,naz+i) = starims(ixx,iyy) if naz-i ge 0 then istar(ixh,iyv,naz-i) = starims(ixx,iyy) ENDFOR ; print, naz-4, ixx,iyy, istar(ixh,iyv,naz-4) ; print, naz-3, ixx,iyy, istar(ixh,iyv,naz-3) ; print, naz-2, ixx,iyy, istar(ixh,iyv,naz-2) ; print, naz-1, ixx,iyy, istar(ixh,iyv,naz-1) ; print, naz, ixx,iyy, istar(ixh,iyv,naz) ; print, naz+1, ixx,iyy, istar(ixh,iyv,naz+1) ; print, naz+2, ixx,iyy, istar(ixh,iyv,naz+2) ; print, naz+3, ixx,iyy, istar(ixh,iyv,naz+3) ; print, naz+4, ixx,iyy, istar(ixh,iyv,naz+4) ; print, naz+5, ixx,iyy, istar(ixh,iyv,naz+5) ixh = ixh + 1 ENDFOR iyv = iyv + 1 ENDFOR naz = naz - iazi ; ; Clean up last bit ; if naz lt 0 then begin iyv = 0 FOR iyy=ivv,ivp do BEGIN ixh = 0 FOR ixx=ihh,ihp do BEGIN FOR i=1,iazd2 do BEGIN if naz+i ge 0 then istar(ixh,iyv,naz+i) = starims(ixx,iyy) ENDFOR ixh = ixh + 1 ENDFOR iyv = iyv + 1 ENDFOR endif ; ENDFOR ; 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 within the boundaries of the windowed ; area ; iset = 0 ; iset = 1 if iset eq 1 then begin ; Do not use filter i1 = ar goto, end1 endif 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-iave-1 gt ivv then ivv = miny ; ivp = nv - iave - 1 if maxy+iave+1 lt ivp then ivp = maxy ; ihh = iave-1 if minx-iave-1 gt ihh then ihh = minx ; ihp = nh - iave - 1 if maxx+iave+1 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=1,iave do BEGIN iix = fix((i*ax)+0.5) iiy = fix((i*ay)+0.5) ixm = ixx - iix iym = iyy - iiy ixp = ixx + iix iyp = iyy + iiy ; print, azb1,azb2,dazbd,azb1,azb2,ixx,iyy,ixp,iyp,ax,ay 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' end1: 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 allor0, bests0,bestam0,bestan0,az1,az2,ht1,ht2,xc,yc,scx,scy, $ minx,maxx,miny,maxy,spvn,ampvn,angvn ; azb1 = az1 htb1 = ht1 azb2 = az2 htb2 = ht2 ; ivv = miny ivp = maxy ihh = minx ihp = maxx ; 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 spvn(ixx,iyy) = bests0 ampvn(ixx,iyy) = bestam0 angvn(ixx,iyy) = bestan0 endif endif ENDFOR ENDFOR 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,expt,imgg,a ; ; Image reader ; this function reads in a LASCO difference image and plots it ; a = float(readfits(name,hdr))/expt ; if imgg eq 1 then begin mxdif = mxbv - mnbv mxd3 = mxdif/3.0 mx = mxbv + mxd3 mn = mnbv - mxd3 print, 'mxbv mnbv mx mn', mxbv, mnbv, mx, mn ; print, ih, iv, !d.name TWIN, /new, xsize=ih+20, ysize=iv 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. ; TWIN, /new, xsize=ih, ysize=iv ; 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)) TWIN, /new, xsize=ih, ysize=iv ; 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)) TWIN, /new, xsize=ih, ysize=iv ; 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,expt,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)/expt ; TWIN, /new, xsize=ih, ysize=iv ; 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)) TWIN, /new, xsize=ih+20, ysize=iv ; 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,expt,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)/expt ; 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)) TWIN, /new, xsize=ih, ysize=iv ; 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 image11,name,ih,iv,az1,az2,ht1,ht2,xc,yc,scx,scy,expt,imgg, $ minx,maxx,miny,maxy,value,badvalue,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)/expt ; 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 ; badcolor = Black 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)) TWIN, /new, xsize=ih, ysize=iv ; 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 va eq badvalue then begin vlu = [badcolor] endif else begin 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 endelse TV, vlu, ixx, iyy endif endif ENDFOR ENDFOR ; endif ; END ; function imagesp, sps,dif,difv,dif0,fract1,fract2,ht1,nspmax,naz, $ psps,pamps,pangs,aznaz1,aznaz2,inis,inif, $ xsizi,ysizi,fspeed,az1,az2,imgg,ifstspflg,iwinnsp,npspnq ; ; This function plots out a contoured speed analysis of speed versus azimuth. ; The units countoured are the amplitudes of the correlation function. ; if imgg eq 1 or imgg eq 3 then begin ; ; MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; print, 'Contouring speed amplitudes' if ifstspflg eq 0 then begin TWIN, /new, xsize=xsizi+80, ysize=ysizi+90 iwinnsp = !d.window - 32 endif if ifstspflg eq 1 then begin TWIN, /show, winnr = iwinnsp endif xsizii = xsizi + 80.0 ysizii = ysizi + 90.0 ifstspflg = 1 ; xticks = (az1-az2)/10.0 xticksi = fix(xticks) xbegin = az1-90.0 if xbegin gt 0.0 then xbegin = fix(xbegin + 0.5) if xbegin lt 0.0 then xbegin = fix(xbegin - 0.5) yticks = fspeed/100.0 yticksi = fix(yticks) ; mini1 = min(sps(*,*)) maxi1 = max(sps(*,*)) print, 'speed image mini1 maxi1', mini1, maxi1 ; printf, one, '1st image mini1 maxi1', mini1, maxi1 sps(*,*) = sps(*,*) - mini1 sps(*,*) = sps(*,*)*(0.5)/(maxi1-mini1) contour, sps(*,*), 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, $ xtickn=strcompress(xbegin-indgen(xticksi+1)*10,/rem), xticks=xticks, $ ; xtickn=['80','70','60','50','40'], xticks=4, $ ; ytickn=['0','100','200','300'], yticks=3, $ ytickn=strcompress(indgen(yticksi+1)*100,/rem), yticks=yticks, $ ; ytickn=['0','100','200','300','400','500','600','700','800','900', $ ; '1000','1100','1200','1300','1400','1500'], yticks=15, $ xstyle=1, ystyle=1, $ pos=[60.0/xsizii,70.0/ysizii,(xsizii-20.0)/xsizii,(ysizii-20.0)/ysizii],$ xtitle='AZIMUTH',ytitle='SPEED' ; mini1 = mini1*100.0 maxi1 = maxi1*100.0 if mini1 ge 0.0 then min = fix(mini1) if mini1 lt 0.0 then min = - fix(abs(mini1)) if maxi1 ge 0.0 then max = fix(maxi1) if maxi1 lt 0.0 then max = - fix(abs(maxi1)) xs = 20 x1 = xsizii-90 yt = ysizii-15 y1 = 25 y2 = 10 ht = fix(ht1) htr = fix((ht1 - ht)*100.0 + 0.5) xyouts, xs, yt, 'SPEED SUMMARY', /device if htr gt 9 then xyouts, [xsizii-120.0], yt, 'HEIGHT = .' , /device if htr le 9 then xyouts, [xsizii-120.0], yt, 'HEIGHT = .0' , /device xyouts, [xsizii-120.0 + 28] , yt, ht , /device xyouts, [xsizii-120.0 + 50] , yt, htr , /device xyouts, x1-22, y1, 'max X 100 =' , /device xyouts, x1+20, y1, max , /device xyouts, x1-22, y2, 'min X 100 =' , /device xyouts, x1+20, y2, min , /device xyouts, 5, yt, strcompress(npspnq,/rem), /device difn = fltarr(naz+1) dif0n = fltarr(naz+1) frac1 = fltarr(naz+1) frac2 = fltarr(naz+1) pspi = fltarr(2) azi = fltarr(2) dbot0 = 0.0 difvn = fltarr(naz+1) dmeana = fltarr(naz+1) dmeana(*) = 0.0 dift = 0.0 diff = 0.0 ndift = 0.0 ndiff = 0.0 for i=0,naz do begin if difv(i) ne 0.0 then begin dift = dift + difv(i) ndift = ndift + 1.0 endif if dif(i) ne 0.0 then begin diff = diff + dif(i) ndiff = ndiff + 1.0 endif endfor if ndift ne 0.0 then difm = dift/ndift if ndiff ne 0.0 then diffm = diff/ndiff difvn = dmeana for i=0,naz do begin if dif0(i) ne 0.0 then begin dif0n(i) = nspmax*50.0*dif0(i)/(fspeed*3.0) + dmeana(i) endif if dif0(i) ne 0.0 and dif(i) ne 0.0 then begin difn(i) = nspmax*50.0*(dif0(i) + dif(i))/(fspeed*3.0) + dbot0 endif if difv(i) ne 0.0 then begin difvn(i) = nspmax*50.0*(difv(i)-difm)/(fspeed*1.0) + dmeana(i) endif frac1(i) = fract1(i)*nspmax*5.0 frac2(i) = fract2(i)*nspmax*5.0 endfor oplot, frac1, linestyle=1 ; fraction of maximum image 1 - dotted oplot, frac2, linestyle=2 ; fraction of maximum image 2 - dashed oplot, dif0n, /noclip ; Difference from monthly minimum / 20.0 oplot, difn, /noclip ; Total diff. from monthly minimum / 20.0 oplot, difvn, /noclip, linestyle=5 ; Difference image - long dashes oplot, dmeana, linestyle=1 ; Mean value - dotted for i=inis, inif do begin azi(0) = (az1 - (aznaz1(i) + aznaz2(i))/2.0) azi(1) = azi(0) + 1.0 pspi(0) = nspmax*psps(i)/fspeed pspi(1) = pspi(0) oplot, azi, pspi endfor endif ; END ; function imageg, angs,dif,difv,dif0,fract1,fract2,ht1,nspmax,naz, $ psps,pamps,pangs,aznaz1,aznaz2,inis,inif, $ xsizi,ysizi,fspeed,az1,az2,imgg,ifstagflg,iwinnag, $ angl,npspnq ; ; This function plots out a contoured angle analysis of angle versus azimuth. ; The units countoured are the angles of the maximum correlation function. ; if imgg eq 1 or imgg eq 3 then begin ; ; MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; print, 'Contouring angles' if ifstagflg eq 0 then begin TWIN, /new, xsize=xsizi+80, ysize=ysizi+90 iwinnag = !d.window - 32 endif if ifstagflg eq 1 then begin TWIN, /show, winnr = iwinnag endif xsizii = xsizi + 80.0 ysizii = ysizi + 90.0 ; xticks = (az1-az2)/10.0 xticksi = fix(xticks) xbegin = az1-90.0 if xbegin gt 0.0 then xbegin = fix(xbegin + 0.5) if xbegin lt 0.0 then xbegin = fix(xbegin - 0.5) yticks = fspeed/100.0 yticksi = fix(yticks) ; ifstagflg = 1 ; mini1 = min(angs(*,*)) mini1 = -angl ; maxi1 = max(angs(*,*)) maxi1 = angl print, 'speed image mini1 maxi1', mini1, maxi1 ; printf, one, '1st image mini1 maxi1', mini1, maxi1 angs(*,*) = angs(*,*) - mini1 angs(*,*) = angs(*,*)*(0.55)/(maxi1-mini1) contour, angs(*,*), levels = [0.0,0.05,0.1,0.15,0.2,0.25,0.3,0.35, $ 0.4,0.45,0.5,0.55], $ c_label = [1,1,1,1,1,1,1,1,1,1,1,1,1], $ c_colors = [121,128,136,144,152,160,120,112,104,96,88,81], /fill, $ xtickn=strcompress(xbegin-indgen(xticksi+1)*10,/rem), xticks=xticks, $ ; xtickn=['80','70','60','50','40'], xticks=4, $ ; ytickn=['0','100','200','300'], yticks=3, $ ytickn=strcompress(indgen(yticksi+1)*100,/rem), yticks=yticks, $ ; ytickn=['0','100','200','300','400','500','600','700','800','900', $ ; '1000','1100','1200','1300','1400','1500'], yticks=15, $ xstyle=1, ystyle=1, $ pos=[60.0/xsizii,70.0/ysizii,(xsizii-20.0)/xsizii,(ysizii-20.0)/ysizii],$ xtitle='AZIMUTH',ytitle='ANGLE' ; if mini1 ge 0.0 then min = fix(mini1) if mini1 lt 0.0 then min = - fix(abs(mini1)) if maxi1 ge 0.0 then max = fix(maxi1) if maxi1 lt 0.0 then max = - fix(abs(maxi1)) xs = 20 x1 = xsizii-90 yt = ysizii-15 y1 = 25 y2 = 10 ht = fix(ht1) htr = fix((ht1 - ht)*100.0 + 0.5) xyouts, xs, yt, 'ANGLE SUMMARY', /device if htr gt 9 then xyouts, [xsizii-120.0], yt, 'HEIGHT = .' , /device if htr le 9 then xyouts, [xsizii-120.0], yt, 'HEIGHT = .0' , /device xyouts, [xsizii-120.0 + 28] , yt, ht , /device xyouts, [xsizii-120.0 + 50] , yt, htr , /device xyouts, x1, y1, 'max =' , /device xyouts, x1+20, y1, max , /device xyouts, x1, y2, 'min =' , /device xyouts, x1+20, y2, min , /device xyouts, 5, yt, strcompress(npspnq,/rem), /device ; iskip = 1 if iskip eq 1 then goto, end1 ; Skip around the overplots on angle dif0n = fltarr(naz+1) frac1 = fltarr(naz+1) frac2 = fltarr(naz+1) dbot0 = 0.0 difvn = fltarr(naz+1) dmeana = fltarr(naz+1) dmeana(*) = 0.0 dift = 0.0 ndift = 0.0 for i=0,naz do begin if difv(i) ne 0.0 then begin dift = dift + difv(i) ndift = ndift + 1.0 endif endfor if ndift ne 0.0 then difm = dift/ndift difvn = dmeana for i=0,naz do begin if difv(i) ne 0.0 then begin difvn(i) = (difv(i)-difm)/1.0 + dmeana(i) endif dif0n(i) = (dif0(i) + difv(i) - difm)/3.0 + dbot0 frac1(i) = fract1(i)*nspmax*5.0 frac2(i) = fract2(i)*nspmax*5.0 endfor oplot, frac1, linestyle=1 ; fraction of maximum image 1 - dotted oplot, frac2, linestyle=2 ; fraction of maximum image 2 - dashed oplot, dif0n, /noclip ; Difference from monthly minimum / 10.0 oplot, difvn, /noclip, linestyle=5 ; Difference image - long dashes oplot, dmeana, linestyle=1 ; Mean value - dotted end1: 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 imagecor, i11s,i22s,ccfs,xxiss,yyjss,xsiz,ysiz,xsizi,ysizi,az1,az2, $ ht1,ht2,hts,spee,amp,ang,dtime,iazv,naz,amaxs,imgg,one, $ ifstflg,iwinn1,iwinn2,iwinn3,iokflg ; ; This function plots the two best images correlated and the correlation. ; The function overplots the specifics of each image on the image and the ; correlation. ; The function plots an "OK" and a "NOT OK" box on the correlation and the ; observer chooses to accept or not accept the correlation. ; if imgg eq 6 or imgg eq 1 then begin if amaxs(iazv) ne -1.0 then begin print, 'set of images correlated and correlation', amaxs(iazv) printf, one, 'set of images correlated and correlation', amaxs(iazv) ; MakeRGB RgbTop = 255 SpaceColors, RgbTop, $ rgbo, reds, greens, blues, $ White, Black, Purple, Red, Green, $ Greyish, Yellowish, Blueish, Reddish, Greenish ; print, 'Contouring first high maximum image' if ifstflg eq 0 then begin TWIN, /new, xsize=xsizi, ysize=ysizi iwinn1 = !d.window - 32 endif if ifstflg eq 1 then begin TWIN, /show, winnr = iwinn1 endif mini1 = min(i11s(*,*,iazv)) maxi1 = max(i11s(*,*,iazv)) print, '1st image mini1 maxi1', mini1, maxi1 ; printf, one, '1st image mini1 maxi1', mini1, maxi1 i11s(*,*,iazv) = i11s(*,*,iazv) - mini1 i11s(*,*,iazv) = i11s(*,*,iazv)*(0.5)/(maxi1-mini1) contour, i11s(*,*,iazv), 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] ; if mini1 ge 0.0 then min = fix(mini1) if mini1 lt 0.0 then min = - fix(abs(mini1)) if maxi1 ge 0.0 then max = fix(maxi1) if maxi1 lt 0.0 then max = - fix(abs(maxi1)) x1 = xsizi-90 yt = ysizi-15 y1 = 25 y2 = 10 az = fix((az1+az2)/2.0 + 0.49) xyouts, [xsizi/4.0], yt, 'AZIMUTH =' , /device xyouts, [xsizi/4.0 + 35] , yt, az , /device ht = fix(ht1) htr = fix((ht1 - ht)*100.0 + 0.5) if htr gt 9 then xyouts, [xsizi/2.0], yt, 'HEIGHT = .' , /device if htr le 9 then xyouts, [xsizi/2.0], yt, 'HEIGHT = .0' , /device xyouts, [xsizi/2.0 + 28] , yt, ht , /device xyouts, [xsizi/2.0 + 50] , yt, htr , /device xyouts, x1, y1, 'max =' , /device xyouts, x1+20, y1, max , /device xyouts, x1, y2, 'min =' , /device xyouts, x1+20, y2, min , /device ; print, 'Contouring second high maximum image' if ifstflg eq 0 then begin TWIN, /new, xsize=xsizi, ysize=ysizi iwinn2 = !d.window - 32 endif if ifstflg eq 1 then begin TWIN, /show, winnr = iwinn2 endif mini2 = min(i22s(*,*,iazv)) maxi2 = max(i22s(*,*,iazv)) print,'2nd image mini2 maxi2', mini2, maxi2 ; printf, one, '2nd image mini2 maxi2', mini2, maxi2 i22s(*,*,iazv) = i22s(*,*,iazv) - mini2 i22s(*,*,iazv) = i22s(*,*,iazv)*(0.5)/(maxi2-mini2) contour, i22s(*,*,iazv), 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] ; if mini2 ge 0.0 then min = fix(mini2) if mini2 lt 0.0 then min = - fix(abs(mini2)) if maxi2 ge 0.0 then max = fix(maxi2) if maxi2 lt 0.0 then max = - fix(abs(maxi2)) xyouts, [xsizi/4.0], yt, 'AZIMUTH =' , /device xyouts, [xsizi/4.0 + 35] , yt, az , /device ht = fix(ht2) htr = fix((ht2 - ht)*100.0 + 0.5) if htr gt 9 then xyouts, [xsizi/2.0], yt, 'HEIGHT = .' , /device if htr le 9 then xyouts, [xsizi/2.0], yt, 'HEIGHT = .0' , /device xyouts, [xsizi/2.0 + 28] , yt, ht , /device xyouts, [xsizi/2.0 + 28] , yt, ht , /device xyouts, [xsizi/2.0 + 50] , yt, htr , /device xyouts, x1, y1, 'max =' , /device xyouts, x1+20, y1, max , /device xyouts, x1, y2, 'min =' , /device xyouts, x1+20, y2, min , /device ; xsizc = xsizi*0.45 ysizc = ysizi*0.45 ; print, 'Contouring correlation of two previous images - OK?' if ifstflg eq 0 then begin TWIN, /new, xsize=xsizc, ysize=ysizc iwinn3 = !d.window - 32 endif if ifstflg eq 1 then begin TWIN, /show, winnr = iwinn3 endif minic = min(ccfs(*,*,iazv)) maxic = max(ccfs(*,*,iazv)) print,'Correlation minic maxic', minic, maxic ; printf, one, 'Correlation minic maxic', minic, maxic cc = fltarr(xsiz,ysiz,naz+1) cc(*,*,iazv) = ccfs(*,*,iazv) - minic cc(*,*,iazv) = cc(*,*,iazv)*(0.5)/(maxic-minic) contour, cc(*,*,iazv), 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] ; if minic ge 0.0 then min = fix(minic*1000.) if minic lt 0.0 then min = fix(abs(minic*1000.)) if maxic ge 0.0 then max = fix(maxic*1000.) if maxic lt 0.0 then max = fix(abs(maxic*1000.)) yt = ysizc-15 xs = 5 spd = fix(spee + 0.5) xyouts, xs, yt, 'sp=', /device xyouts, [xs+1], yt, spd, /device xam = 65 ap = fix(amp*1000. + 0.5) xyouts, xam, yt, 'amp=0.', /device xyouts, [xam+15], yt, ap, /device xan = 135 if ang ge 0.0 then an = fix(ang + 0.5) if ang lt 0.0 then an = -fix(abs(ang) + 0.5) xyouts, xan, yt, 'ang=', /device xyouts, [xan+3], yt, an, /device ht = fix(hts) htr = fix((hts - ht)*100.0 + 0.5) if htr eq 100 then begin ht = ht + 1 htr = 0 endif x1 = xsizc-100 ya = ysizc-35 y0 = ysizc-50 y1 = ysizc-65 y2 = ysizc-80 dt = dtime/60.0 dti = fix(dt) dtr = fix((dt - dti)*100.0 + 0.5) if dtr eq 100 then begin dt = dt + 1 dtr = 0 endif if dtr gt 9 then xyouts, x1, ya, ' dt = .' , /device if dtr le 9 then xyouts, x1, ya, ' dt = .0' , /device if dti ge 10 then xyouts, [x1 + 10] , ya, dti , /device if dti lt 10 then xyouts, [x1 + 10] , ya, dti , /device xyouts, [x1 + 28] , ya, dtr , /device if htr gt 9 then xyouts, x1, y0, ' ht = .' , /device if htr le 9 then xyouts, x1, y0, ' ht = .0' , /device xyouts, [x1 + 10] , y0, ht , /device xyouts, [x1 + 28] , y0, htr , /device if maxic ge 0.0 then xyouts, x1, y1, 'max = 0.' , /device if maxic lt 0.0 then xyouts, x1, y1, 'max = -0.' , /device xyouts, x1+30, y1, max , /device if minic ge 0.0 then xyouts, x1, y2, 'min = 0.' , /device if minic lt 0.0 then xyouts, x1, y2, 'min = -0.' , /device xyouts, x1+30, y2, min , /device ; ; ; Plot a cross on the correlation at the location of the maximum ; nmid = (xsiz-1.0)/2.0 xyouts, xxiss(iazv)+nmid, yyjss(iazv)+nmid, '+' ; y4 = 15 xyouts, [xsizc-20], y4, 'OK?' ,/device ; ; Draw two boxes on the bottom right of the correlation image plot ; ; Box 2 is the "NOT OK" box. ; ih = xsizc-20 iv = 35 abox2 = fltarr(5) aboy2 = fltarr(5) abox2(0) = ih+2 abox2(1) = ih+2 abox2(2) = ih+18 abox2(3) = ih+18 abox2(4) = ih+2 aboy2(0) = 2 + iv aboy2(1) = 18 + iv aboy2(2) = 18 + iv aboy2(3) = 2 + iv aboy2(4) = 2 + iv ; ; Box 1 is the "OK" box. ; iv = 0 abox1 = fltarr(5) aboy1 = fltarr(5) abox1(0) = ih+2 abox1(1) = ih+2 abox1(2) = ih+18 abox1(3) = ih+18 abox1(4) = ih+2 aboy1(0) = 2 + iv aboy1(1) = 18 + iv aboy1(2) = 18 + iv aboy1(3) = 2 + iv aboy1(4) = 2 + iv ; oplot, abox1,aboy1, thick = 2, /noclip oplot, abox2,aboy2, thick = 2, /noclip ; ok = 0 iokflg = getok(xsizc-20,ysizc,ok) ; ifstflg = 1 ; endif 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,expt ; ; 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) aa = readfits(name1,hdraa) expt = fxpar(hdraa,'EXPTIME') print, ' Exposure time', expt ; 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,instp=instp,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,iflg=iflg, $ xsize=xsize,ysize=ysize,filtriave=iave,nimage=nimage,nsmo=nsmo, $ astar=astar,iline=iline,iedge=iedge,msts=msts,mstl=mstl,mste=mste, $ mmi=mmi ; ; 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 or a predefined set of heights and azimuth angles, the ; researcher defines an area on the image from which to determine ; outward motion. This area is drawn on the next image and is ; also defined with a default outward motion on this image. The LASCO ; data is filtered using the best available filter, and the following ; image within that same azimuth in the data set is, too. 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 at a number ; of predetermined heights before continuing. The best answer within ; given 1 deg azimuth areas are listed and plotted as the best ; correct answer. The structures are NOT tracked. The program redefines ; the area first defined on the first image and speeds, amplitudes and angles ; within these areas for each pixel are drawn, coded according to the ; highest amplitude values of the different speeds between sspeed and fspeed ; maped to the pixels flagged. At the end of processing two images in the ; set, the program uses the second of the images in the set as the first, ; and then uses the next image in the set as the pair and processes these ; two images in a similar manner over the height specified. At the end of ; the processing at a given height, the program can repeat at a new height ; if the user so specifies. ; 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 imgst = img if n_elements(mmi) eq 0 then mmi = 0 ; month min img pre = 1 if n_elements(istep) eq 0 then istep = 4 ; Image pair seperation isteps = istep if n_elements(instp) eq 0 then instp = istep ; 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 value ; if n_elements(dhts) eq 0 then dhts = 1.0 ; Height step multiplier if n_elements(az1) eq 0 then az1 = 35 ; bigger azimuth (1) if n_elements(az2) eq 0 then az2 = 25 ; lesser azimuth (2) if n_elements(ht1) eq 0 then ht1 = 8.0 ; smaller height (1) if n_elements(ht2) eq 0 then ht2 = 8.5 ; 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 astar = 0 ; no star correlation iline = 0 ; no line correlation iedge = 0 ; no edge correlation if n_elements(msts) ne 0 then astar = 1 ; a star correlation if n_elements(mstl) ne 0 then iline = 1 ; a line correlation if n_elements(mste) ne 0 then iedge = 1 ; an edge correlation if astar eq 0 then msts = 1.0 ; define msts if iline eq 0 then mstl = 1.0 ; define mstl if iedge eq 0 then mste = 1.0 ; define mste 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 if n_elements(iflg) eq 0 then iflg = 1 ; if 0, no corr. max plt 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 itrigger2 = 0 ; dir1 = '\LASCO Data\C3 Data\' dir2 = '\LASCO Data\C3 Results\' fini1 = '.dif' fini2 = '.dat' fini3 = '.sum' openw, one, dir1+name+strcompress(img,/rem)+fini2,/get_lun openw, two, dir2+name+strcompress(img,/rem)+fini3,/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 = 4000 ini = 0 ; istar = fltarr(2*xsize+1,2*ysize+1,naz+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 npspnq = 1 imageseq = 0 ; if nimage eq 3 then imageseq = 5 if nimage eq 3 then imageseq = 3 ; npspnq = 4 ; if nimage eq 4 then imageseq = 7 if nimage eq 4 then imageseq = 4 ; if nimage eq 5 then imageseq = 9 if nimage gt 5 then imageseq = 5 if imageseq eq 0 then imageseq = nimage pcos = 1.0 psin = 0.0 pang = 0.0 ht1s = ht1 ht2s = ht2 htdn = dhts*(ht2s - ht1s) azdn = 0.0 az1s = az1 az2s = az2 psp = pspd ifstflg = 0 ifstspflg = 0 ifstagflg = 0 firsttimetru = 0 inis=0 ; 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,expt1) ; 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 = 200/expt1 mxbv1 = 40/expt1 mnbv1 = 20/expt1 ; mnbv1 = 160/expt1 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 ; ; Mar 09 if nam1 eq 'mar09_' then begin mxbv1 = 40/expt1 mnbv1 = -10/expt1 hcent1 = 518.784 vcent1 = 18.989 endif ; up: img1 = image(namea,hsize1,vsize1,mxbv1,mnbv1,expt1,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) if firsttimetru eq 0 then begin az1naz = fltarr(naz+1) az2naz = fltarr(naz+1) bestsp0 = fltarr(naz+1) bestamp0 = fltarr(naz+1) bestang0 = fltarr(naz+1) ndifv = fltarr(naz+1) difv = fltarr(naz+1) dif = fltarr(naz+1) dif0 = fltarr(naz+1) diffs = fltarr(naz+1) dif00s = fltarr(naz+1) difvs = fltarr(naz+1) diffss = fltarr(naz+1) dif00ss = fltarr(naz+1) difvss = fltarr(naz+1) ndifvs = 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 endif ; ; 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 ht1 ht2',az1,az2,ht1/rs,ht2/rs 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,expt2) 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/expt2 mnbv2 = 20/expt2 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 ; ; Mar 09 if nam1 eq 'mar09_' then begin mxbv2 = 40/expt2 mnbv2 = -10/expt2 hcent2 = 518.784 vcent2 = 18.989 endif ; fd=findmx(az11,az22,ht11,ht22,hcent2,vcent2,scx2,scy2,xsize,ysize,iave, $ minx11,maxx11,miny11,maxy11) print, minx11,maxx11,miny11,maxy11 ; ; Plot the second differenced image on a window with the new area ; superimposed. ; img2 = image(nameb,hsize2,vsize2,mxbv2,mnbv2,expt2,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 ; ; If a difference minus the monthly minimum image is available for use, ; read it in. The name of this image is assumed in the same directory ; as the data, and it is named the same as the difference data, but the ; img value is 0 and the step is also 0. Thus, the January 15 monthly ; minimum data file would be named jan15_0.dif ; if firsttimetru eq 0 then begin mm0 = fltarr(hsize1,vsize1) if mmi eq 1 then begin imgq = 0 istepq = 0 par = param(name,imgq,istepq,fini1,name0,hsize0,vsize0,hcent0,vcent0,$ scx0,scy0,day0,time0,expt0) mxmn = mxmnbv(b,hsize0,vsize0,az11,az22,ht11,ht22,hcent0,vcent0, $ scx0,scy0,minx11,maxx11,miny11,maxy11,mxbv0,mnbv0) mm0 = fltarr(hsize0,vsize0) ; mxbv0 = 40.0/expt0 ; mnbv0 = -10.0/expt0 img0 = image(name0,hsize0,vsize0,mxbv0,mnbv0,expt0,imgg,mm0) endif endif ; ; 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) ; ; Provide a difference image from the base smoothed and with stars removed ; for plotting later. ; if mmi eq 1 then mm = (ar + br)/2.0 ; if firsttimetru eq 0 then begin istar = fltarr(2*xsize+1,2*ysize+1,naz+1) endif istar(*,*,*) = -10000.0 ; if astar eq 0 then goto, nostar ; ; Pseudo stellar images ; if firsttimetru eq 0 then begin stari = fltarr(hsize1,vsize1) endif stari(*,*) = 1.0 stput = starput(stari,hsize1,vsize1,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2) ; smo = smoo(stari,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,starir) if nsmo ge 3 then begin ; ; If not special C2 filtered data, remove next line ; starir = starir - smooth(starir,15) starir = smooth(starir,nsmo) print, 'Stellar images smoothed',nsmo,' by',nsmo printf, one, 'Stellar images smoothed',nsmo,' by',nsmo endif ; ft1 = filtr(starir,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,starims) ; ; Extract stellar images from data to use as patterns to correlate with ; ; if firsttimetru eq 0 then begin ; starims = fltarr(hsize1,vsize1) ; endif ; starims(*,*) = starir(*,*) ; exim = extrims(starims,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,istar) ; nostar: ; if firsttimetru eq 0 then begin ilnh = fltarr(2*xsize+1,2*ysize+1,naz+1) ilnv = fltarr(2*xsize+1,2*ysize+1,naz+1) endif ilnh(*,*,*) = -10000.0 ilnv(*,*,*) = -10000.0 ; if iline eq 0 then goto, noline ; ; Pseudo horizontal line images ; if firsttimetru eq 0 then begin ilnri = fltarr(hsize1,vsize1) endif ilnri(*,*) = 1.0 lnput = hlineput(ilnri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2) ; smo = smoo(ilnri,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,ilnrir) if nsmo ge 3 then begin ilnrir = smooth(ilnrir,nsmo) print, 'Horizontal line images smoothed',nsmo,' by',nsmo printf, one, 'Horizontal line images smoothed',nsmo,' by',nsmo endif ; ft1 = filtr(ilnrir,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,linims) ; ; Extract horizontal line images from data to use as patterns to correlate ; exim = extrims(linims,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,ilnh) ; ; Pseudo vertical line images ; if firsttimetru eq 0 then begin ilnri = fltarr(hsize1,vsize1) endif ilnri(*,*) = 1.0 lnput = vlineput(ilnri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2) ; smo = smoo(ilnri,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,ilnrir) if nsmo ge 3 then begin ilnrir = smooth(ilnrir,nsmo) print, 'Vertical line images smoothed',nsmo,' by',nsmo printf, one, 'Vertical line images smoothed',nsmo,' by',nsmo endif ; ft1 = filtr(ilnrir,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,linims) ; ; Extract vertical line images from data to use as patterns to correlate ; exim = extrims(linims,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,ilnv) ; noline: ; if firsttimetru eq 0 then begin ieh = fltarr(2*xsize+1,2*ysize+1,naz+1) iev = fltarr(2*xsize+1,2*ysize+1,naz+1) endif ieh(*,*,*) = -10000.0 iev(*,*,*) = -10000.0 ; if iedge eq 0 then goto, noedge ; ; Pseudo horizontal edge images ; if firsttimetru eq 0 then begin iedgri = fltarr(hsize1,vsize1) endif iedgri(*,*) = 1.0 edgput = hedgeput(iedgri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2) ; smo = smoo(iedgri,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,iedgrir) if nsmo ge 3 then begin iedgrir = smooth(iedgrir,nsmo) print, 'Horizontal edge images smoothed',nsmo,' by',nsmo printf, one, 'Horizontal edge images smoothed',nsmo,' by',nsmo endif ; ft1 = filtr(iedgrir,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,edgims) ; ; Extract horizontal edge images from data to use as patterns to correlate ; exim = extrims(edgims,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,ieh) ; ; Pseudo vertical edge images ; if firsttimetru eq 0 then begin iedgri = fltarr(hsize1,vsize1) endif iedgri(*,*) = 1.0 edgput = vedgeput(iedgri,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2) ; smo = smoo(iedgri,hsize1,vsize1,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,iedgrir) if nsmo ge 3 then begin iedgrir = smooth(iedgrir,nsmo) print, 'Vertical edge images smoothed',nsmo,' by',nsmo printf, one, 'Vertical edge images smoothed',nsmo,' by',nsmo endif ; ft1 = filtr(iedgrir,hsize1,vsize1,iave,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,minx,maxx,miny,maxy,edgims) ; ; Extract vertical edge images from data to use as patterns to correlate ; exim = extrims(edgims,xsize,ysize,hcent1,vcent1,scx1,scy1, $ az1,az2,ht1,ht2,iev) ; noedge: ; ; 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 if firsttimetru eq 0 then begin spv = fltarr(hsize1,vsize1) ampv = fltarr(hsize1,vsize1) angv = fltarr(hsize1,vsize1) spvn = fltarr(hsize1,vsize1) ampvn = fltarr(hsize1,vsize1) angvn = fltarr(hsize1,vsize1) xxis = fltarr(hsize1,vsize1) yyjs = fltarr(hsize1,vsize1) endif spv(*,*) = badspv ampv(*,*) = badamp angv(*,*) = badang spvn(*,*) = badspv ampvn(*,*) = badamp angvn(*,*) = badang if firsttimetru eq 0 then begin ccfst = fltarr(xsize,ysize,nk+1) nccfst = fltarr(hsize1,vsize1) amaxb = fltarr(hsize1,vsize1) endif goodcorr = 1.0 amaxb(*,*) = goodcorr if firsttimetru eq 0 then begin i11s = fltarr(2*xsize+1,2*ysize+1,naz+1) i22s = fltarr(2*xsize+1,2*ysize+1,naz+1) ccfs = fltarr(xsize,ysize,naz+1) htbb1 = fltarr(naz+1) htbb2 = fltarr(naz+1) xxiss = fltarr(naz+1) yyjss = fltarr(naz+1) endif spcw = (xsize-1.0)*scx1/dtime if npspnq eq 1 then begin ; sperpix1 = spcw/(xsize-1.0) nper100 = fix(100.0/sperpix1) if nper100 gt 0 then begin spbin = 100.0/nper100 endif else begin spbin = 200.0 endelse nspmax = fix(fspeed/spbin) nspmin = fix(sspeed/spbin) if firsttimetru eq 0 then begin spsums = fltarr(naz+1,nspmax+1) fract1 = fltarr(naz+1) fract2 = fltarr(naz+1) endif spsums(*,*) = 0.0 diffs(*) = 0.0 dif00s(*) = 0.0 difvs(*) = 0.0 ndifvs(*) = 0.0 fract1(*) = 0.0 fract2(*) = 0.0 if firsttimetru eq 0 then begin badsum = fltarr(naz+1) endif badsum(*) = 0.0 if firsttimetru eq 0 then begin angsums = fltarr(naz+1,nspmax+1) endif angsums(*,*) = 0.0 if firsttimetru eq 0 then begin nsums = fltarr(naz+1,nspmax+1) endif nsums(*,*) = 0.0 if firsttimetru eq 0 then begin spsumss = fltarr(nspmax+1) endif spsumss(*) = 0.0 if firsttimetru eq 0 then begin angsumss = fltarr(nspmax+1) endif angsumss(*) = 0.0 if firsttimetru eq 0 then begin nsumss = fltarr(nspmax+1) endif nsumss(*) = 0.0 endif if ht1 ne ht1s then begin spsums(*,*) = 0.0 angsums(*,*) = 0.0 nsums(*,*) = 0.0 spsumss(*) = 0.0 angsumss(*) = 0.0 nsumss(*) = 0.0 endif ; if firsttimetru eq 0 then begin spsum = fltarr(nspmax+1) endif spsum(*) = 0.0 if firsttimetru eq 0 then begin spsu = fltarr(naz+1,nspmax+1) endif spsu(*,*) = 0.0 if firsttimetru eq 0 then begin sps = fltarr(naz+1,nspmax+1) endif sps(*,*) = 0.0 if firsttimetru eq 0 then begin angsum = fltarr(nspmax+1) endif angsum(*) = 0.0 if firsttimetru eq 0 then begin angsu = fltarr(naz+1,nspmax+1) endif angsu(*,*) = 0.0 if firsttimetru eq 0 then begin angs = fltarr(naz+1,nspmax+1) endif angs(*,*) = 0.0 if firsttimetru eq 0 then begin nsum = fltarr(nspmax+1) endif nsum(*) = 0.0 if firsttimetru eq 0 then begin spsumsss = fltarr(nspmax+1) angsumsss = fltarr(nspmax+1) endif spsumsss(*) = 0.0 angsumsss(*) = 0.0 diffs(*) = 0.0 difvs(*) = 0.0 dif00s(*) = 0.0 diffss(*) = 0.0 dif00ss(*) = 0.0 difvss(*) = 0.0 ndifvs(*) = 0.0 if firsttimetru eq 0 then begin amaxs = fltarr(naz+1) amaxss = fltarr(naz+1) bestht = fltarr(naz+1) endif amaxs(*) = -1.0 amaxss(*) = -1.0 ; 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 bestamp0(*) = 0.0 bestang0(*) = 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). ; amaxb(*,*) = goodcorr badsum(*) = 0.0 sxsize3 = spcw/3.0 ; 1/3 the speed equivalent to the length of xsize psp = sspeed dif(*) = 0.0 dif0(*) = 0.0 difv(*) = 0.0 ndifv(*) = 0.0 again: print, ' psp', psp printf, one, ' psp', psp ; ccfs(0,0,*) = -10000.0 amaxss(*) = amaxs(*) xxis(*,*) = -999.0 yyjs(*,*) = -999.0 ; angln = angl sspeedn = sspeed fspeedn = fspeed if psp eq 0.0 then begin angln = 360.0 sspeedn = 0.0 ; fspeedn = 40.0 endif 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,amaxss,i11s,i22s,ccfs,istar,ilnh,ilnv, $ ieh,iev,ccfst,nccfst,amaxb,mggg,badsum, $ az11,az22,ht11,ht22,htbb1,htbb2,xxiss,yyjss,xsize,one, $ angln,sspeedn,fspeedn,smax,badspv,badang,badamp, $ spv,angv,ampv,xxis,yyjs,msts,mstl,mste) ; ; Determine the highest amplitude correlations of speed and angle values ; within the area defined by az1naz, az2naz and ht1, ht2. ; iiii = 0 spsumiilim = 1.0 for i=0,naz do begin if psp eq 0.0 and badsum(i) eq 0.0 then goto, skip if psp ne 0.0 and sspeed eq 0.0 then begin if badsum(i) eq 0.0 then goto, skip if nsums(i,0) eq 0.0 then goto, skipa1 if spsums(i,0)/nsums(i,0) gt spsumiilim then goto, skip skipa1: endif bests0 = bestsp0(i) bestam0 = bestamp0(i) bestan0 = bestang0(i) spsum(*) = 0.0 angsum(*) = 0.0 nsum(*) = 0.0 ; fd = findmx(az1naz(i),az2naz(i),ht1,ht2,hcent1,vcent1,scx1,scy1, $ xsize,ysize,iave,mnx,mxx,mny,mxy) ; ; Contour the highest amplitude correlations per azimuth angle and check ; each to see if it is a valid correlation. If so flag it as good and ; allow it as a best value for this correlation set. ; xsizi = 450.0 ysizi = xsizi iazv = i iokflg = 0 ; Correlation bad unless found otherwise or imager disabled ; ; The following analysis gives the best values within the windowed azimuth and ; height area triggered by the best valid amplitude value within that area ; bv=bestvn(spv,ampv,angv,hsize1,vsize1,hcent1,vcent1,scx1,scy1,xsize,ysize, $ az1naz(i),az2naz(i),ht1,ht2,mnx,mxx,mny,mxy,bests0,bestam0,bestan0) ; ; The following compiles an amplitude sum keyed on velocity over the range of ; velocity from sspeed to fspeed within the windowed azimuth and height area. ; The sum uses all valid values within the windowed area. ; sv=sumv(spv,ampv,angv,i1,i2,mm,mm0,ccfst,nccfst,spcw,spbin,xxis,yyjs, $ hsize1,vsize1,hcent1,vcent1,iazv,scx1,scy1,xsize,ysize, $ az1naz(i),az2naz(i),iazv,ht1,ht2,mnx,mxx,mny,mxy,spsum,angsum, $ diff,difvv,dif00,ndifvv,nsum,nspmax,nspmin,iflg,psp, $ angl,sspeed,fspeed,smax) ; if mmi eq 1 then begin if psp gt 0.0 then begin if ndifvv gt 0.0 then begin dif(i) = diff difv(i) = difvv dif0(i) = dif00 ndifv(i) = ndifvv diffs(i) = diffs(i) + dif(i) dif00s(i) = dif00s(i) + dif0(i) difvs(i) = difvs(i) + difv(i) ndifvs(i) = ndifvs(i) + ndifv(i) ; ; Standard deviation values of the images for the best correlated images ; xsized2 = xsize/2 txsizemd2 = xsize*2-xsized2 ysized2 = ysize/2 tysizemd2 = ysize*2-ysized2 mi11s = 0.0 mi22s = 0.0 ni12s = 0.0 FOR ix=xsized2,txsizemd2 DO BEGIN FOR iy=ysized2,tysizemd2 DO BEGIN mi11s = i11s(ix,iy,i) + mi11s mi22s = i22s(ix,iy,i) + mi22s ni12s = ni12s + 1.0 ENDFOR ENDFOR mi11s = mi11s/ni12s mi22s = mi22s/ni12s ; difi1mean = 0.0 difi2mean = 0.0 FOR ix=xsized2,txsizemd2 DO BEGIN FOR iy=ysized2,tysizemd2 DO BEGIN difi1mean = (i11s(ix,iy,i)-mi11s)*(i11s(ix,iy,i)-mi11s) +difi1mean difi2mean = (i22s(ix,iy,i)-mi22s)*(i22s(ix,iy,i)-mi22s) +difi2mean ENDFOR ENDFOR if (ni12s-1.0) gt 0.0 then stdev1 = sqrt(difi1mean/(ni12s-1.0)) if (ni12s-1.0) gt 0.0 then stdev2 = sqrt(difi2mean/(ni12s-1.0)) fract11 = ndifvs(i)*stdev1/dif00s(i) if fract11 lt 0.2 then fract1(i) = fract11 fract22 = ndifvs(i)*stdev2/dif00s(i) if fract22 lt 0.2 then fract2(i) = fract22 ; print, dif00s(i)/ndifvs(i),stdev1,stdev2,fract1(i),fract2(i) ; printf, one, dif00s(i)/ndifvs(i),stdev1,stdev2,fract1(i),fract2(i) endif endif endif ; FOR ii=0,nspmax do BEGIN if psp eq 0.0 and badsum(i) eq 0.0 then goto, skip if psp ne 0.0 and sspeed eq 0.0 then begin if badsum(i) eq 0.0 then goto, skip if nsums(i,0) eq 0.0 then goto, skipa2 if spsums(i,0)/nsums(i,0) gt spsumiilim then goto, skip skipa2: endif iiii = iiii + 1 if nsum(ii) eq 0.0 then goto, skipa3 spsums(i,ii) = spsums(i,ii) + spsum(ii) angsums(i,ii) = angsums(i,ii) + angsum(ii) nsums(i,ii) = nsums(i,ii) + nsum(ii) spsumss(ii) = spsumss(ii) + spsum(ii) angsumss(ii) = angsumss(ii) + angsum(ii) nsumss(ii) = nsumss(ii) + nsum(ii) skipa3: ENDFOR ; ; print, 'nspmax', nspmax ; print, 'spsums', spsums(i,*) ; print, 'angsums', angsums(i,*) ; print, 'nsums', nsums(i,*) ; htl = htbb1(iazv)/rs htu = htbb2(iazv)/rs hts = htl + (bests0*dtime)/rs bestht(i) = hts ; if amaxss(i) eq amaxs(i) and bests0 eq bestsp0(i) then begin print, format='(A9,I4,4F8.2,8X,2F8.2,F7.3)', $ ' bad ccf ', i,az1naz(i),az2naz(i),ht1/rs,ht2/rs, $ bestsp0(i),bestang0(i),bestamp0(i) endif ; if amaxss(i) ne amaxs(i) and bests0 ne bestsp0(i) then begin print, format='(A9,I4,7F8.2,F7.3)', $ ' good ccf', i,az1naz(i),az2naz(i),ht1/rs,ht2/rs,hts, $ bests0,bestan0,bestam0 ; iokflg = 1 ; ; Plot out the best correlation of the bunch ; ; noplot = 0 noplot = 1 if noplot ne 0 then begin cr=imagecor(i11s,i22s,ccfs,xxiss,yyjss,xsize,ysize,xsizi,ysizi, $ az1naz(i),az2naz(i),htl,htu,hts,bests0,bestam0,bestan0,dtime, $ iazv,naz,amaxss,imgg,one,ifstflg,iwinn1,iwinn2,iwinn3,iokflg) ; ; noplot = 0 ; if noplot ne 0 then begin if astar ne 0 then begin cr=imagecor(i11s,istar,ccfs,xxiss,yyjss,xsize,ysize,xsizi,ysizi, $ az1naz(i),az2naz(i),htl,htu,hts,bests0,bestam0,bestan0,dtime, $ iazv,naz,amaxss,imgg,one,ifstflg,iwinn1,iwinn2,iwinn3,iokflgs) endif ; if iline ne 0 then begin cr=imagecor(i11s,ilnh,ccfs,xxiss,yyjss,xsize,ysize,xsizi,ysizi, $ az1naz(i),az2naz(i),htl,htu,hts,bests0,bestam0,bestan0,dtime, $ iazv,naz,amaxss,imgg,one,ifstflg,iwinn1,iwinn2,iwinn3,iokflgl) ; cr=imagecor(i11s,ilnv,ccfs,xxiss,yyjss,xsize,ysize,xsizi,ysizi, $ az1naz(i),az2naz(i),htl,htu,hts,bests0,bestam0,bestan0,dtime, $ iazv,naz,amaxss,imgg,one,ifstflg,iwinn1,iwinn2,iwinn3,iokflgl) endif ; if iedge ne 0 then begin cr=imagecor(i11s,ieh,ccfs,xxiss,yyjss,xsize,ysize,xsizi,ysizi, $ az1naz(i),az2naz(i),htl,htu,hts,bests0,bestam0,bestan0,dtime, $ iazv,naz,amaxss,imgg,one,ifstflg,iwinn1,iwinn2,iwinn3,iokflge) ; cr=imagecor(i11s,iev,ccfs,xxiss,yyjss,xsize,ysize,xsizi,ysizi, $ az1naz(i),az2naz(i),htl,htu,hts,bests0,bestam0,bestan0,dtime, $ iazv,naz,amaxss,imgg,one,ifstflg,iwinn1,iwinn2,iwinn3,iokflge) endif endif ; endif ; if not iokflg then goto, skip ; if iokflg set to 0 then skip next ; iokflg = 0 ; iokflag reset to bad ; ; print, az1naz(i),az2naz(i),ht1/rs,ht2/rs,bestam0,bests0,bestan0 ; ; The following is supposed to allow only the good maxima through ; amaxs(i) = bestam0 ; bestsp0(i) = bests0 bestamp0(i) = bestam0 bestang0(i) = bestan0 ; ; Turn arrays to plot into all the same or nothing values at a given azimuth ; alon = allor0(bests0,bestam0,bestan0,az1naz(i),az2naz(i),ht1,ht2, $ hcent1,vcent1,scx1,scy1,mnx,mxx,mny,mxy,spvn,ampvn,angvn) ; skip: endfor ; if psp gt fspeed then goto, next if iiii eq 0 then goto, nextimg 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, 'amplitude0 =', bestamp0 printf, one, 'amplitude0 =', bestamp0 ; print, 'angle0 =', bestang0 printf, one, 'angle0 =', bestang0 ; ; Determine the mean value of the speeds and determine the standard ; deviation of the sample ; avesp = 0.0 aveht = 0.0 nspds = 0.0 for i=0,naz do begin avesp = avesp + bestsp0(i) if bestsp0(i) gt 30.0 then aveht = aveht + bestht(i) if bestsp0(i) gt 30.0 then nspds = nspds + 1.0 endfor if nspds ne 0.0 then avesp = avesp/nspds if nspds ne 0.0 then aveht = aveht/nspds difmean2 = 0.0 for i=0,naz do begin if bestsp0(i) gt 30.0 then begin difmean2 = (bestsp0(i) - avesp)*(bestsp0(i) - avesp) + difmean2 endif endfor stdev = 0.0 devmean = 0.0 if (nspds-1.0) gt 0.0 then stdev = sqrt(difmean2/(nspds-1.0)) if nspds gt 1.0 then devmean = stdev/sqrt(nspds) print, 'image # speeds mean speed std sdevmean ave height' print, format='(f5.0,f7.0,2x,4f12.2)', img, nspds, avesp, stdev, devmean, aveht printf, one, 'image # speeds mean speed std sdevmean ave height' printf, one, format='(f5.0,f7.0,2x,4f12.2)', img, nspds, avesp, stdev, devmean, aveht if firsttimetru eq 0 then $ printf, two, 'image # speeds mean speed std sdevmean ave height' printf, two, format='(f5.0,f7.0,2x,4f12.2)', img, nspds, avesp, stdev, devmean, aveht firsttimetru = 1 ; ; 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 ; img11 = image11(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ expt,imgg,minx,maxx,miny,maxy,spvn,badspv,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 ; img11 = image11(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ expt,imgg,minx,maxx,miny,maxy,ampvn,badamp,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 ; img11 = image11(namea,hsize1,vsize1,az1,az2,ht1,ht2,hcent1,vcent1,scx1,scy1, $ expt,imgg,minx,maxx,miny,maxy,angvn,badang,-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) ; ; Plot the speed and angle sums ; spsu(*,*) = 0.0 angsu(*,*) = 0.0 spsumsss(*) = 0.0 angsumsss(*) = 0.0 FOR i=0,naz do BEGIN FOR ii=0,nspmax do BEGIN if nsums(i,ii) ne 0.0 then begin spsu(i,ii) = spsums(i,ii)/nsums(i,ii) angsu(i,ii) = angsums(i,ii)/nsums(i,ii) endif ENDFOR print, 'naz nspmax', i, nspmax ; print, 'spsums', spsums(i,*) ; print, 'angsums', angsums(i,*) print, format = '(1x,a4,/,12f6.3,/,12f6.3,/,12f6.3)','spsu',spsu(i,*) print, format = '(1x,a4,/,12f6.1,/,12f6.1,/,12f6.1)','angsu',angsu(i,*) print, format = '(1x,a4,/,12f6.0,/,12f6.0,/,12f6.0)','nsums',nsums(i,*) ENDFOR FOR ii=0,nspmax do BEGIN if nsumss(ii) ne 0.0 then begin spsumsss(ii) = spsumss(ii)/nsumss(ii) angsumsss(ii) = angsumss(ii)/nsumss(ii) endif ENDFOR FOR i=0,naz do BEGIN FOR ii=0,nspmax do BEGIN sps(i,ii) = spsu(i,ii) angs(i,ii) = angsu(i,ii) ENDFOR if ndifvs(i) ne 0.0 then begin diffss(i) = diffs(i)/ndifvs(i) difvss(i) = difvs(i)/ndifvs(i) dif00ss(i) = dif00s(i)/ndifvs(i) endif ENDFOR ; print, 'nspmax ht1 ht2 spsumsss', nspmax, ht1/rs, ht2/rs, spsumsss printf, one, 'nspmax ht1 ht2 spsumsss', nspmax, ht1/rs, ht2/rs, spsumsss printf, two, 'nspmax ht1 ht2 spsumsss', nspmax, ht1/rs, ht2/rs, spsumsss print, 'nspmax ht1 ht2 angsumsss', nspmax, ht1/rs, ht2/rs, angsumsss printf, one, 'nspmax ht1 ht2 angsumsss', nspmax, ht1/rs, ht2/rs, angsumsss printf, two, 'nspmax ht1 ht2 angsumsss', nspmax, ht1/rs, ht2/rs, angsumsss print, 'nspmax ht1 ht2 nsumss', nspmax, ht1/rs, ht2/rs, nsumss printf, one, 'nspmax ht1 ht2 nsumss', nspmax, ht1/rs, ht2/rs, nsumss printf, two, 'nspmax ht1 ht2 nsumss', nspmax, ht1/rs, ht2/rs, nsumss ; ; ; Do the correlations once more if npspnq lt nimage ; for i=0,naz do begin pspi = bestsp0(i) pang = bestang0(i) pamp = bestamp0(i) ; ; determine heights and azimuths to be used as the best guess on this ; set of 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 ; ; 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 ; xsizim = naz*6.0 ysizim = fspeed/3.0 ; inif = ini - 1 ; imgsa = imageg(angs,diffss,difvss,dif00ss,fract1,fract2,htm/rs,nspmax,naz, $ psps,pamps,pangs,aznaz1,aznaz2,inis,inif, $ xsizim,ysizim,fspeed,az1,az2,imgg,ifstagflg,iwinnag,angl, $ npspnq) ; imgs = imagesp(sps,diffss,difvss,dif00ss,fract1,fract2,htm/rs,nspmax,naz, $ psps,pamps,pangs,aznaz1,aznaz2,inis,inif, $ xsizim,ysizim,fspeed,az1,az2,imgg,ifstspflg,iwinnsp,npspnq) ; print, 'npspnq nimage', npspnq, nimage ; nextimg: ; npspnq = npspnq + 1 ; Bump image seq. count ; if npspnq gt imageseq then begin ; Update parameters ; inis = ini az1 = az1 + azdn ; New starting azimuth 1 az2 = az2 + azdn ; New starting azimuth 2 ht1s = ht1 ht1 = ht1 + htdn ; New starting height 1 ht2 = ht2 + htdn ; New starting height 2 ; 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) printf, 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 ; if ht1 gt htmax then begin ; if height beyond limit, then end ; 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) printf, 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 ; goto, end1 ; end correlations ; endif ; istep = isteps npspnq = 1 img = imgst ifstflg = 0 ifstspflg = 0 ifstagflg = 0 ; goto, start ; endif ; ; We are not finished with the image sequence ; if nimage eq 3 then begin if npspnq eq 4 then begin istep = 2*isteps img = imgst - 2*isteps endif if npspnq eq 5 then img = imgst endif ; if nimage eq 4 then begin if npspnq eq 5 then begin istep = 2*isteps img = imgst - 2*isteps endif if npspnq eq 7 then img = imgst endif ; if nimage eq 5 then begin if npspnq eq 6 then begin istep = 2*isteps img = imgst - 2*isteps endif if npspnq eq 8 then img = imgst endif ; instp = istep ; img = img + instp ; New image ; 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) printf, 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 free_lun, two ; end