;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; This procedure helps determine the velocity of solar matter from ;; eruptive prominences. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Here is a little helper function to calculate the distance of particles ;; in solar radii. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION SOLAR, radius, x_cen, y_cen, x_loc, y_loc dist_sq=( (x_loc-x_cen)^2 + (y_loc-y_cen)^2 ) solar_dist= SQRT(dist_sq) / radius; return, solar_dist end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Here is a little helper function to convert radians into degrees. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION RAD_2_DEG, rads IF ( rads LT 0 ) THEN $ degs = ( rads/!PI ) * 180 + 360 $ ELSE $ degs = ( rads/!PI ) * 180 return, degs end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Here is a little helper function to determine the angle measured in ;; degrees ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION CALC_ANGLE, sun_x, sun_y, _X_, _Y_ dX = _X_- sun_x & dY = _Y_- sun_y IF (dX NE 0) THEN BEGIN IF (dY EQ 0) THEN BEGIN IF (dX GT 0) THEN theta_rad = 0.0 $ ELSE theta_rad = !PI ENDIF ELSE $ theta_rad = ATAN(dY, dX) ENDIF ELSE BEGIN $ IF (dY GT 0) THEN theta_rad = !PI / 2.0 $ ELSE theta_rad = -(!PI / 2.0) ENDELSE theta_deg = RAD_2_DEG(theta_rad) return, theta_deg end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Here is a little helper function to determine the position angle in ;; solar coordinates. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; FUNCTION GET_PA, angle_meas, correction_angle if ( correction_angle LT 0 ) then begin if ( angle_meas LT (-correction_angle) ) then $ PA = (360 + correction_angle) + angle_meas $ else $ PA = angle_meas + correction_angle endif else begin if ( (360 - angle_meas ) LT correction_angle ) then $ PA = correction_angle - angle_meas $ else $ PA = angle_meas + correction_angle endelse return, PA end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; This procedure will compute the velocity of the "point" between ;; two successive measurements. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO VEL_CALC, last_position, last_angle, last_time, $ curr_position, curr_angle, curr_time, $ speed, sol_dist if ( last_angle GE curr_angle ) then phi = last_angle - curr_angle $ else phi = curr_angle - last_angle phi = phi * !PI / 180 dist_sq = (last_position^2) + (curr_position^2) - $ 2.0 * last_position * curr_position * cos(phi) sol_dist = sqrt(dist_sq) print,FORMAT='("Distance moved is: ", F)', sol_dist print,FORMAT='("Angle subtended is: ", F)', phi print,FORMAT='("Current time: ", I6, "Last time: ", I6)',curr_time, last_time if (curr_time LT last_time) then curr_time = curr_time + 240000 del_time = curr_time - last_time case 1 of (del_time GT 100) AND (del_time LT 200): del_time = del_time - 100 + 20 (del_time GT 200) AND (del_time LT 300): del_time = del_time - 200 + 40 (del_time GT 300) AND (del_time LT 400): del_time = del_time - 300 + 60 (del_time GT 400) AND (del_time LT 500): del_time = del_time - 400 + 80 (del_time GT 500) AND (del_time LT 600): del_time = del_time - 500 + 100 else: del_time = del_time - 40 endcase print,FORMAT='("Diff. in time is: ", I)', del_time speed = 0.0D speed = (sol_dist * 6.9599 * 10^5) / del_time print,FORMAT='("Velocity is: ", D)', speed return end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; This procedure will print the distance travelled and the velocity of ;; mass in question. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO P_VEL, file_num, solar_distance, velocity printf,file_num, FORMAT='(36X, F, F)', solar_distance, velocity end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; This procedure prints the correct info. to an output file in a ;; specified format. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO P_2_FILE, file_num, sol_dist, pos_ang, time printf,file_num, FORMAT='(F, F, I12)', sol_dist, pos_ang, time end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; This procedure prints the correct info. to an output file in a ;; format useful with the fit program for velocity calculations etc.. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO P_4_FIT, file_num, year, day, time, sol_rad, sol_ang tmh = strmid(time,0,4) tsec = strmid(time,5,2) tmh = fix(tmh) & tsec = fix(tsec) printf,file_num, FORMAT='(I6, I5, I6, I4, F, F3, F, F3, A)', $ year, day, tmh, tsec, sol_rad, 0.05, sol_ang, 0.25, "dPmon" end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Now for the main program...... ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Set up a good color table.... w/b ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; _a=intarr(256) & _b=intarr(256) & _c=intarr(256) openr,15, '/home/puu/cordyn/color_idl/wb0.rgb' readu,15,_a & readu,15,_b & readu,15,_c tvlct,_a,_b,_c close,15 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Install a pretty cool cursor ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; device, cursor_image=[0,1,2,4,40,48,56,64,0,64,56,48,40,4,2,1],cursor_xy=[15,7] pixX=675 & pixY=900 image=bytarr(pixX,pixY) answer=' ' filenm=' ' infile='vel_list' outfile=' ' print,FORMAT='($,"Enter the name of the output file: ")' read, outfile print, ' ' ;; There will be two possible format outputs.... print,FORMAT='($,"Will you be using (fit) format(1), or standard(2): ")' read,style if (stlye = 1) then begin print,FORMAT='($,"Enter the year ")' read,img_year print,FORMAT='($,"Enter the DOY ")' read, img_day img_feature=' ' print,FORMAT='($,"Enter the feature discription ")' read, img_feature endif QUITTING = 0 ;;;;;;;; Variable used for controlling the while loop... ;; First, find out how many images there are.... close,1 & openr,1,infile num=0 while ( not eof(1) ) do begin readf,1,filenm & num=num+1 endwhile close,1 imgnum=num ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; The following structure will be used to hold the information ;; specific to each image for use in printing to an output file in a ;; specified format. The structure simplifies tracking when moving ;; back and forth through the images. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; plt_pnt = { POI, do_print: 0, year: 0, day: 0, hr_min: 0, sec: 0, $ sol_rad: 0.0, sol_rad_err: 0.0, sol_ang: 0.0, $ sol_ang_err: 0.0, TYPE: ' ' } ;; NOW, MAKE AN ARRAY OF THIS STRUCTURE.......... IMG_PLT = REPLICATE( {POI}, imgnum ) IMG_PLT.do_print = 0 IMG_PLT.year = img_year IMG_PLT.day = img_day IMG_PLT.sol_rad_err = 0.05 IMG_PLT.sol_ang_err = 0.25 IMG_PLT.TYPE = 'dPmon' IMG_POINTER = 0 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Using the first image ( and therefore assuming correct alignment ) the ;; solar center will be "determined" (there is margin for error here!!!!) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; close,2 & openr,2,infile & readf,2,filenm close,3 & openr,3,filenm & readu,3,image & close,3 window, 0, xsize=pixX, ysize=pixY, title='Solar Center Test' wset,0 tv, bytscl(image,min=0,max=255) print,'Please click the mouse as close to the solar' print,'center as possible.' print,' ' cursor,_x_temp, _y_temp, /DOWN, /DEV ledge = _x_temp & WHILE(image(ledge,_y_temp) LT 100) $ DO ledge=ledge-1 redge = _x_temp & WHILE(image(redge,_y_temp) LT 100) $ DO redge=redge+1 top = _y_temp & WHILE(image(_x_temp, top) LT 100) $ DO top=top+1 bottom = _y_temp & WHILE(image(_x_temp, bottom) LT 100) $ DO bottom=bottom-1 XDIAM = redge - ledge YDIAM = top - bottom print, FORMAT='("XDIAM is: ", I)', XDIAM print, FORMAT='("YDIAM is: ", I)', YDIAM if ( XDIAM GE YDIAM ) then begin DIAM = XDIAM sol_rad = DIAM/2 sol_cen_x = redge - fix(DIAM/2) sol_cen_y = top - fix(DIAM/2) endif else begin DIAM = YDIAM sol_rad = DIAM/2 sol_cen_x = redge - fix(DIAM/2) sol_cen_y = top - fix(DIAM/2) endelse print, FORMAT='("The Diameter is: ", I)', DIAM print, FORMAT='("The solar radius is: ", I)', sol_rad print, FORMAT='("The center is at x= ", I ," y= ", I)', sol_cen_x, sol_cen_y print, ' ' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; We must also determine the solar angle (South Pole is measured here) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; readf,2,filenm & close,2 close,3 & openr,3,filenm & readu,3,image & close,3 window, 0, xsize=pixX, ysize=pixY, title='Solar South Pole Test' tv, bytscl(image,min=0,max=255) print,'Please click the mouse as close to the solar' print,'south pole as possible. Angle, not distance, is' print,'most important here.' print,' ' cursor,x_south,y_south,/DOWN, /DEV SP_angle = CALC_ANGLE(sol_cen_x, sol_cen_y, x_south, y_south) theta_corr = 180 - SP_angle print,FORMAT='("The South Pole angle is: ", F)', SP_angle print,FORMAT='("This gives a correction angle of: ", F)', theta_corr print,' ' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; Done finding center of sun and the South Pole angle... ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; infile_pos=lonarr(imgnum) && infile_tmp=long(0) & ifp=0 outfile_pos=lonarr(imgnum) && outfile_tmp=long(0) & ofp=0 close,11 & openr,11,infile if (style = 2) then begin close,21 & openw,21,outfile printf,21,FORMAT='(A15,A15,A10)','Solar Radius', 'Position Angle', 'Time' endif print,"Use the left mouse button to advance to next image." print,"Use the center mouse button to mark a point on the image." print,"Use the right mouse button to go back one image." print,' ' num_2_print = 0 while( NOT QUITTING ) do begin point_lun,-11,infile_tmp & infile_pos(ifp)=infile_tmp if (style = 2) then $ point_lun,-21,outfile_tmp & outfile_pos(ofp)=outfile_tmp readf,11,filenm print,FORMAT='("Current file is: ", A6)', filenm close,31 & openr,31,filenm & readu,31,image & close,31 shrt_name = strmid(filenm,4,4) if ( shrt_name EQ '14.a' ) then $ tv,bytscl(image,min=0, max=190) $ else if ( shrt_name EQ '47.a' ) then $ tv,bytscl(image,min=0, max=225,top=240) $ else $ tv,bytscl(image,min=0, max=215,top=240) xyouts, 0.05,0.95, filenm, size=2, /NORMAL cursor,_x,_y, /down,/dev & mouse_button = !err if ( mouse_button EQ 1 ) then begin ; Go Forward.... ifp=ifp+1 & ofp=ofp+1 & IMG_POINTER = IMG_POINTER+1 if ( ifp GT imgnum-1 ) then begin ifp=ifp-1 & ofp=ofp-1 & IMG_POINTER = IMG_POINTER-1 point_lun,11,infile_pos(ifp) point_lun,21,outfile_pos(ofp) print,'-------------- AT END -------------------' print, FORMAT='($,"##### DO YOU WANT TO QUIT(y/n)? " )' read,answer if (answer = 'y') then QUITTING = 1 endif endif else begin if ( mouse_button EQ 2 ) then begin ; Get x, y and advance to next image print,'Click the mouse button at desired location...' print,' ' ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; This little bit of code is used to display the pixel location of the ;; cursor along with the value of the pixel for use in tracking a particular ;; "mass" from image to image. The code is "borrowed" with slight modification ;; from the rdpix.pro routine of idl. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; s = size(image) s(1) = s(1)-1 s(2) = s(2)-1 !err = 0 if s(s(0)+1) GE 4 then form = 'F' else form = 'I' cr = string("15b) form="($,'x=',i4,', y=',i4,', value=',"+form+", a)" while (!err LT 1) do begin tvrdc,x_pos,y_pos,2,/dev if (x_pos LE s(1)) and (y_pos LE s(2)) and (x_pos GE 0) and $ (y_pos GE 0) then begin if (!order EQ 1) then y_pos = s(2) - y_pos else y_pos = y_pos print,form = form, x_pos,y_pos,image(x_pos,y_pos),cr endif endwhile ;; Here ends the "borrowed" code............ print,FORMAT='("The last intensity was: ", I, A)',image(x_pos,y_pos),cr true_dist = SOLAR(sol_rad, sol_cen_x, sol_cen_y, x_pos, y_pos) meas_angle = CALC_ANGLE(sol_cen_x, sol_cen_y, x_pos, y_pos) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; The meas_angle is not in SOLAR COORDINATES, therefore the angle must ;; be adjusted to properly determint the Position Angle. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; position_angle = GET_PA(meas_angle, theta_corr) print,FORMAT='("The PA measured is: ", F)', position_angle print,' ' if (style = 2) then $ print,'Writing to file......' $ else $ print,'Storing data.........' print,' ' if (style = 2) then begin num_2_print = num_2_print + 1 if ( num_2_print GT 1 ) then begin VEL_CALC, last_pos, last_ang, last_time, $ true_dist, position_angle, curr_time, speed, sol_dist P_VEL, 21, sol_dist, speed endif P_2_FILE, 21, true_dist, position_angle, filenm endif else begin THM = strmid(filenm,0,4) & THM = fix(THM) TSEC = strmid(filenm,4,2) & TSEC = fix(TSEC) IMG_PLT(IMG_POINTER).do_print = 1 IMG_PLT(IMG_POINTER).hr_min = THM IMG_PLT(IMG_POINTER).sec = TSEC IMG_PLT(IMG_POINTER).sol_rad = true_dist IMG_PLT(IMG_POINTER).sol_ang = position_angle endelse ifp=ifp+1 & ofp=ofp+1 & IMG_POINTER=IMG_POINTER+1 if ( ifp GT imgnum-1 ) then begin ifp=ifp-1 & ofp=ofp-1 & IMG_POINTER=IMG_POINTER-1 point_lun,11,infile_pos(ifp) if (style = 2) then $ point_lun,21,outfile_pos(ofp) print,' ' print,'------------ AT END / AT END ---------' print,' ' print, FORMAT='($,"##### DO YOU WANT TO QUIT(y/n)? " )' read,answer if (answer = 'y') then QUITTING = 1 endif endif endelse if ( mouse_button EQ 4 ) then begin ; GO BACKWARD print,'------------Going Backward---------------' ifp = ifp-1 & if (ifp LE 0) then begin ifp=0 & print, 'At Beginning....' & print,' ' endif point_lun,11,infile_pos(ifp) ofp = ofp-1 & if (ofp LE 0) then ofp=0 if (style = 2) then $ point_lun,21,outfile_pos(ofp) IMG_POINTER = IMG_POINTER - 1 & if (IMG_POINTER LT 0) then $ IMG_POINTER = 0 endif endwhile close,11 if (style = 2) then begin printf,21,FORMAT='( I )', 1 printf,21,"y" printf,21,"y" printf,21,"n" close,21 endif else begin close,21 & openw,21,outfile printf, 21, outfile printf, 21, FORMAT='(F4.2,F6.2)', 1.0, 2.5 printf, 21, img_feature data_points = 0 for num_of_points = 0, imgnum-1 do begin if ( IMG_PLT(num_of_points).do_print ) then data_points = data_points + 1 endfor printf, 21, FORMAT='(I2)', data_points for count = 0, imgnum-1 do begin if ( IMG_PLT(count).do_print ) then $ printf, 21, FORMAT='(I4,I5,I6,I4,4F8.2, A6)', IMG_PLT(count).year, $ IMG_PLT(count).day, IMG_PLT(count).hr_min, $ IMG_PLT(count).sec, IMG_PLT(count).sol_rad, $ IMG_PLT(count).sol_rad_err, IMG_PLT(count).sol_ang, $ IMG_PLT(count).sol_ang_err, IMG_PLT(count).TYPE endfor printf, 21, FORMAT='(I2)', 1 printf, 21,"y" printf, 21,"y" printf, 21,"n" endelse end