wcs_demo.pro

basic median filter simulation

  printf,file_unit,"if (poles[0] ne -1) then $"
  printf,file_unit,$
    "  latitude_inv(poles) = latitude_inv(poles)/abs(latitude_inv(poles))*9.d1"
endif

printf,file_unit, $
    "dist = sphdist(longitude,latitude,longitude_inv,latitude_inv,/degrees)"
printf,file_unit,"erase"
printf,file_unit,$
"print,'The largest arrow on the plot will represent a difference of '"
printf,file_unit,"print,max(dist),' degrees.'"
printf,file_unit,"read,'Press return to continue',key"
printf,file_unit,$
       "norm = sqrt((longitude-longitude_inv)^2 + (latitude-latitude_inv)^2)"
printf,file_unit,"lon_diff=dist*(longitude-longitude_inv)"
printf,file_unit,"good = where(norm ne 0.d0)"
printf,file_unit,"lon_diff(good) = lon_diff(good)/norm(good)"
printf,file_unit,"lat_diff = dist*(latitude-latitude_inv)"
printf,file_unit,"lat_diff(good) = lat_diff(good)/norm(good)"
printf,file_unit,"velovect,lon_diff,lat_diff,longitude[*,0],latitude[0,*]"
end

; PROCEDURE FOR OPTION 3
pro wcs_rot,file_unit,map,param1,param2
printf,file_unit,";PLOTTING"
printf,file_unit,"; Plot the resulting map."
if ((map ge 0) and (map le 22)) then begin
  ; For all but the spherical cube projections, simply plot the results from
  ; wcssph2xy.pro as is.
  printf,file_unit,"xdelta = (max(xx) - min(xx))/20"
  printf,file_unit,"ydelta = (max(y) - min(y))/20"
  printf,file_unit,$
         "plot,xx,y,psym = 3,xrange = [min(xx) - xdelta,max(xx) + xdelta],$"
  printf,file_unit,$
         "yrange = [min(y) - ydelta,max(y) + ydelta],xstyle = 4,ystyle = 4"
  printf,file_unit,"zero_ind = where(latitude[0,*] eq min(abs(latitude[0,*])))"
  printf,file_unit,$
  "xyouts,xx(lon_index,zero_ind[0]),y(lon_index,zero_ind[0]),$"
  printf,file_unit,$
  "       strcompress(string(long(longitude(lon_index,zero_ind[0])))),$"
  printf,file_unit,"       alignment = 0.5"
  printf,file_unit,$
  "zero_ind2 = where(longitude[*,0] eq min(abs(longitude[*,0])))"
  printf,file_unit,$
  "xyouts,xx(zero_ind2[0],lat_index),y(zero_ind2[0],lat_index),$"
  printf,file_unit,$
  "       strcompress(string(long(latitude(zero_ind2[0],lat_index)))),$"
  printf,file_unit,"       alignment = 0.5"
  printf,file_unit,$
  "non_zero_ind = where(longitude[*,0] ne min(abs(longitude[*,0])))
  printf,file_unit,$
  "for i = 0,zero_ind[0] - 1 do $"
  printf,file_unit,$
  "    oplot,xx(non_zero_ind,i),y(non_zero_ind,i),psym=4"
  printf,file_unit,$
  "for i = zero_ind[0] + 1,n_elements(longitude[0,*]) - 1 do $"
  printf,file_unit,"   oplot,xx(non_zero_ind,i),y(non_zero_ind,i),psym=4"
endif else begin
  printf,file_unit,"xx = -x"
  printf,file_unit,"yy = y"

  printf,file_unit,""
  printf,file_unit,"; Make arrays with the locations of all points."
  printf,file_unit,"face_0 = where(face eq 0)"
  printf,file_unit,"face_1 = where(face eq 1)"
  printf,file_unit,"face_2 = where(face eq 2)"
  printf,file_unit,"face_3 = where(face eq 3)"
  printf,file_unit,"face_4 = where(face eq 4)"
  printf,file_unit,"face_5 = where(face eq 5)"

  printf,file_unit,""
  printf,file_unit,"; Define the size of the box around each face."
  if (map eq 23) then begin
    printf,file_unit,"x_len = 90"
    printf,file_unit,"y_len = 90"
  endif else begin
    printf,file_unit,"x_len = 2*!radeg"
    printf,file_unit,"y_len = 2*!radeg"
  endelse

  printf,file_unit,""
  printf,file_unit,$
  "; Correctly adjust the x and y values for display purposes (they all start "
  printf,file_unit,$
  "; out on the same face)."
  printf,file_unit,"if (face_0[0] ne -1) then begin"
  printf,file_unit,"  x0 = -x(face_0)"
  printf,file_unit,"  y0 = y(face_0) - y_len"
  printf,file_unit,"  xx(face_0) = x0"
  printf,file_unit,"  yy(face_0) = y0"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_1[0] ne -1) then begin"
  printf,file_unit,"  x1 = -x(face_1) + 2.d0*x_len"
  printf,file_unit,"  y1 = y(face_1)"
  printf,file_unit,"  xx(face_1) = x1"
  printf,file_unit,"  yy(face_1) = y1"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_2[0] ne -1) then begin"
  printf,file_unit,"  x2 = -x(face_2) + x_len"
  printf,file_unit,"  y2 = y(face_2)"
  printf,file_unit,"  xx(face_2) = x2"
  printf,file_unit,"  yy(face_2) = y2"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_3[0] ne -1) then begin"
  printf,file_unit,"  x3 = -x(face_3)"
  printf,file_unit,"  y3 = y(face_3)"
  printf,file_unit,"  xx(face_3) = x3"
  printf,file_unit,"  yy(face_3) = y3"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_4[0] ne -1) then begin"
  printf,file_unit,"  x4 = -x(face_4) - x_len"
  printf,file_unit,"  y4 = y(face_4)"
  printf,file_unit,"  xx(face_4) = x4"
  printf,file_unit,"  yy(face_4) = y4"
  printf,file_unit,"endif"
  printf,file_unit,"if (face_5[0] ne -1) then begin"
  printf,file_unit,"  x5 = -x(face_5)"
  printf,file_unit,"  y5 = y(face_5) - y_len"
  printf,file_unit,"  xx(face_5) = x5"
  printf,file_unit,"  yy(face_5) = y5"
  printf,file_unit,"endif"

  printf,file_unit,""
  printf,file_unit,$
  "; Define plot ranges by finding which faces are actually used."
  printf,file_unit,"if (face_4[0] ne -1) then x_low = -3*x_len/2 $"
  printf,file_unit,"else if (face_3[0] ne -1) then x_low = -x_len/2 $"
  printf,file_unit,"else if (face_2[0] ne -1) then x_low = x_len/2 $"
  printf,file_unit,$
  "else if ((face_1[0] ne -1) or (face_5[0] ne -1) or (face_0[0] ne -1)) $"
  printf,file_unit,"  then x_low = 3*x_len/2"
  printf,file_unit,"if (face_4[0] ne -1) then x_high = -x_len/2 $"
  printf,file_unit,"else if (face_2[0] ne -1) then x_high = 3*x_len/2 $"
  printf,file_unit,"else if (face_3[0] ne -1) then x_high = x_len/2 $"
  printf,file_unit,$
  "else if ((face_1[0] ne -1) or (face_5[0] ne -1) or (face_0[0] ne -1)) $"
  printf,file_unit,"  then x_high = 5*x_len/2"
  printf,file_unit,"if (face_5[0] ne -1) then y_low = -3*y_len/2 $"
  printf,file_unit,$
  "else if ((face_4[0] ne -1) or (face_3[0] ne -1) or (face_2[0] ne -1) or $"
  printf,file_unit,"   (face_1[0] ne -1)) then y_low = -y_len/2 $"
  printf,file_unit,"else if (face_0[0] ne -1) then y_low = y_len/2"
  printf,file_unit,"if (face_0[0] ne -1) then y_high = 3*y_len/2 $"
  printf,file_unit,$
  "else if ((face_1[0] ne -1) or (face_3[0] ne -1) or (face_2[0] ne -1) or $"
  printf,file_unit,"   (face_4[0] ne -1)) then y_high = y_len/2 $"
  printf,file_unit,"else if (face_5[0] ne -1) then y_high = -y_len/2"

  printf,file_unit,""
  printf,file_unit,"; Plot the points calculated by wcssph2xy."
  printf,file_unit,$
  "plot,xx,yy,psym=3,xrange=[x_low,x_high],yrange=[y_low,y_high],xstyle=4,$"
  printf,file_unit,"     ystyle=4"
  printf,file_unit,"zero_ind = where(latitude[0,*] eq min(abs(latitude[0,*])))"
  printf,file_unit,$
  "xyouts,xx(lon_index,zero_ind[0]),yy(lon_index,zero_ind[0]),$"
  printf,file_unit,$
  "       strcompress(string(long(longitude(lon_index,zero_ind[0])))),$"
  printf,file_unit,"       alignment = 0.5"
  printf,file_unit,$
  "zero_ind2 = where(longitude[*,0] eq min(abs(longitude[*,0])))"
  printf,file_unit,$
  "xyouts,xx(zero_ind2[0],lat_index),yy(zero_ind2[0],lat_index),$"
  printf,file_unit,$
  "       strcompress(string(long(latitude(zero_ind2[0],lat_index)))),$"
  printf,file_unit,"       alignment = 0.5"
  printf,file_unit,$
  "non_zero_ind = where(longitude[*,0] ne min(abs(longitude[*,0])))
  printf,file_unit,$
  "for i = 0,zero_ind[0] - 1 do $"
  printf,file_unit,$
  "    oplot,xx(non_zero_ind,i),yy(non_zero_ind,i),psym=4"
  printf,file_unit,$
  "for i = zero_ind[0] + 1,n_elements(longitude[0,*]) - 1 do $"
  printf,file_unit,"   oplot,xx(non_zero_ind,i),yy(non_zero_ind,i),psym=4"
endelse
end

; MAIN DEMO PROGRAM
pro wcs_demo
print,""
print,"This demo program demonstrates the basic usage of the IDL procedures"
print,"wcssph2xy.pro and wcsxy2sph.pro.  You will be prompted for information"
print,"about the type of map projection you would like to try out and what"
print,"portion of the sky you would like to view.  All of the commands"
print,"actually issued to carry out these operations will be recorded in a"
print,"journal file so that the user may later reproduce the results from this"
print,"demo by issuing the commands him/herself.  Enjoy!"
key=''
print,""
repeat read,"Enter 'c' to continue or 'x' to exit:",key $
until ((key eq 'c') or (key eq 'x'))
if (key eq 'x') then stop
print,""

; Major loop of whole program.
repeat begin

print,""
print,"Your options are:"
print,"(1) Convert spherical (sky) coordinates to x and y coordinates"
print,"    (in other words, perform a map projection) and plot the results."
print,"(2) Do (1) without plotting, then perform the inverse operation."
print,"    Plot the results, then plot the difference between the original"
print,"    sky coordinates and the coordinates that have been produced by"
print,"    wcssph2xy and wcsxy2sph."
print,"(3) Do (1) with an added twist, rotating the coordinate system."
print,"(4) Exit"
print,""
repeat read,"Enter a number between 1 and 4:",option $
until ((option ge 1) and (option le 4))
print,""

if (option eq 4) then stop

file_name = ""
repeat begin
  read,"Please enter a name for the journal file:",file_name
  print,""
  suffix = strmid(file_name,strlen(file_name)-4,4)
  if (suffix ne ".pro") then file_name = string(file_name,".pro")
  file_test = file_search(file_name)
  if (file_test[0] ne "") then begin
    print,"The file ",file_name," already exists."
    print,"You can overwrite this file, but if you used this journal name"
    print,"previously in this IDL session, you will not get the desired"
    print,"results.  To avoid any conflicts, either quit and start a new"
    print,"session of IDL using this name (and ignore this message) or give a"
    print,"new name to the journal file.  NOTE:  This is due to IDL's"
    print,"inability to re-compile a procedure except from the interactive"
    print,"mode." 
    print,""
    read,"Type 'y' to overwrite the file:",key
    if (key ne 'y') then file_name = ""
  endif
endrep until (file_name ne "")
openw,file_unit,file_name,/get_lun

printf,file_unit,$
"; This is an IDL procedure created by running the IDL program wcs_demo.pro"
printf,file_unit,$
"; and can be executed from the IDL prompt by typing .run ",file_name,"."
printf,file_unit,$
"; This procedure may be far more complicated than what you need.  In order"
printf,file_unit,$
"; to make it more user-friendly, I have broken up the tasks performed into"
printf,file_unit,"; the following categories:"
printf,file_unit,";   (1) SET-UP -- sections declaring constants"
printf,file_unit,$
";   (2) CONVERSION -- section in which spherical to xy conversion is done"
printf,file_unit,$
";   (3) LABELS -- sections setting up and printing labels on the maps"
printf,file_unit,$
";   (4) PLOTTING -- sections in which data or lines are plotted"
printf,file_unit,$
";To find the appropriate section, simply search for one of these four"
printf,file_unit,";capitalized words."

printf,file_unit,""
printf,file_unit,string("pro ",strmid(file_name,0,strlen(file_name) - 4))

map = 0
print,""
print,"Which map projection would you like to try?  Your options are:"
print,"Number  Description                Number  Description"
print,"------  -------------------------  ------  -------------------------"
print,"   0    Default = Cartesian           1    Zenithal perspective"
print,"   2    Gnomic                        3    Orthographic"
print,"   4    Stereographic                 5    Zenithal Equidistant"
print,"   6    Zenithal polynomial (not implemented)"
print,"   7    Zenithal equal area           8    Airy"
print,"   9    Cylindrical perspective      10    Cartesian"
print,"  11    Mercator                     12    Cylindrical equal area"
print,"  13    Conical perspective          14    Conical equidistant"
print,"  15    Conical equal area           16    Conical orthomorphic"
print,"  17    Bonne's equal area           18    Polyconic"
print,"  19    Sanson-Flmsteed              20    Parabolic"
print,"  21    Hammer-Aitoff                22    Mollweide"
print,"  23    Cobe Quadrilateralized Spherical Cube"
print,"  24    Quadrilateralized Spherical Cube"
print,"  25    Tangential Spherical Cube"
print,""
print,$
"NOTE: This demo program does not support the map types: 1-4,8-9,11,13, or 16 "
print,$
"with coordinate system rotation (option 3 above).  These are allowed by"
print,$
"wcssph2xy.pro and wcsxy2sph.pro, but due to problems with the general case of"
print,$
"latitude and longitude restrictions, these map types were skipped here." 
print,""
repeat read,"Please enter a number from 0 to 25:",map $
until ((map ge 0) and (map le 25))

if (option eq 3) then begin
  if ((map le 4) or (map eq 8) or (map eq 9) or (map eq 11) or (map eq 13) $
      or (map eq 16)) then begin
    close,file_unit
    file_delete, file_name
    message,"The map type selected is not supported with coordinate rotations."