wcs_demo.pro

basic median filter simulation

  endif else begin
    print,$
    "The idea behind the rotation of the coordinate systems is to relocate the"
    print,$
  "'special' point of the projection.  For instance, the azimuthal projections"
    print,$
    "project from the north pole.  So, the lines of longitude appear as rays"
    print,$
    "coming from the center of the projection and lines of latitude appear as"
    print,$
    "concentric rings around the center.  By rotating the coordinate system,"
    print,$
    "a different point can play the role of the north pole in this example."
    print,$
    "To perform the rotation, the latitude and longitude of the new 'special'"
    print,$
    "point must be given.  In addition, to specify a full rotation, a third"
    print,$
    "angle must be given.  This angle specifies the longitude of the north"
    print,$
    "pole in the transformed system and has a default of 180 degrees."
    print,""
    read,"Please enter the longitude of the 'special' point:",alpha 
    read,"Please enter the latitude of the 'special' point:",delta
    read,"Please enter the third angle (enter 180 for the default):",longpole
  endelse
endif
  
printf,file_unit,";SET-UP"
printf,file_unit,"; Set-up constants used later in this procedure"
printf,file_unit,"map = ",map
print,""

; Get parameters for map types that require them.
case map of
  1:begin
    read,$
    "AZP: Enter distance of source to projection (range = [0,10^14]):",param1
  end
  6:begin
    close,file_unit
    file_delete,file_name,/allow
    message,"ZPN: This map projection has not been implemented."
  end
  8:begin
    print,"AIR: Enter the angular distance from the tangent point in which the"
    read,"error is to be minimized (range = [0,90]):",param1
  end
  9:begin
    read,"CYP: Enter the radius of the cylinder (range = [0,10^14]):",param2
    print,"CYP: Enter the distance from the projection point to the center of"
    read,"the sphere (range = [-10^14,10^14], but not -radius):",param1
  end
  12:begin
    print,"CEA: Enter the square of the cosine of the latitude at which the"
    read,"map is conformal (range = [0,1]):",param1
  end
  13:begin
    read,$
    "COP: Lower angle at which cone intersects sphere (range = [-90,upper]):",$
    theta1
    read,$
    "COP: Upper angle at which cone intersects sphere (range = [lower,90]):",$
    theta2
    param1 = (theta2+theta1)/2.
    param2 = abs(theta2 - theta1)/2
  end
  14:begin
    read,$
    "COD: Lower angle at which cone intersects sphere (range = [-90,upper]):",$
    param1
    read,$
    "COD: Upper angle at which cone intersects sphere (range = [lower,90]):",$
    param2
  end
  15:begin
    read,$
    "COE: Lower angle at which cone intersects sphere (range = [-90,upper]):",$
    param1
    read,$
    "COE: Upper angle at which cone intersects sphere (range = [lower,90]):",$
    param2
  end
  16:begin
    read,$
    "COO: Lower angle at which cone intersects sphere (range = [-90,upper]):",$
    param1
    read,$
    "COO: Upper angle at which cone intersects sphere (range = [lower,90]):",$
    param2
  end
  17:begin
    read,"BON: Characteristic angle (range = [-90,90]):",param1
  end
  else:
endcase

if (n_elements(param1) ne 0) then printf,file_unit,"param1 = ",param1
if (n_elements(param2) ne 0) then printf,file_unit,"param2 = ",param2
if (n_elements(alpha) ne 0) then printf,file_unit,"alpha = ",alpha
if (n_elements(delta) ne 0) then printf,file_unit,"delta = ",delta
if (n_elements(longpole) ne 0) then printf,file_unit,"longpole = ",longpole

print,"Would you like to:"
print,"(1) Do a whole-sky map."
print,"(2) Select a (rectangular) region on the sky to map."
print,""
repeat read,"Enter '1' or '2':",choice until ((choice eq 1) or (choice eq 2))
print,""

; Set-up to do a full-sky map.
if (choice eq 1) then begin
  ; set-up the longitude range
  printf,file_unit,"min_lon = 0"
  printf,file_unit,"max_lon = 345"
  printf,file_unit,"lon_spacing = 15"

  ; set-up the latitude range (this differs from map to map because some maps
  ; diverge at particular latitudes)
  if ((map eq  1) or (map eq 3)) then begin
    printf,file_unit,"min_lat = 0"
    printf,file_unit,"max_lat = 90"
  endif else if (map eq 2) then begin
    printf,file_unit,"min_lat = 15"
    printf,file_unit,"max_lat = 90"
  endif else if (map eq 4) then begin
    printf,file_unit,"min_lat = -75"
    printf,file_unit,"max_lat = 90"
  endif else if (map eq 8) then begin
  ; For the Airy map projection, the minimum usable latitude depends on the
  ; input parameters, so it must be calculated now.
    xi = (findgen(90) + 1)/!radeg
    xi_b = (!pi/2.0 - param1/!radeg)/2.0
    radius=-!radeg*(alog(cos(xi))/tan(xi)+alog(cos(xi_b))/tan(xi_b)*tan(xi))
    i = 0
    repeat i = i + 1 $
    until ((radius[i + 1] lt radius[i]) or (i eq (n_elements(radius) - 2)))
    i = i - 1
    min_lat = 90 - 2*!radeg*xi[i]
    printf,file_unit,"min_lat = ",min_lat[0]
    printf,file_unit,"max_lat = 90"
  endif else if (map eq 9) then begin
 ; The CYP map projection diverges at the poles when param1 (mu) is equal to 0.
    if (param1 eq 0) then begin
      printf,file_unit,"min_lat = -75"
      printf,file_unit,"max_lat = 75"
    endif else begin
      printf,file_unit,"min_lat = -90"
      printf,file_unit,"max_lat = 90"
    endelse
  endif else if (map eq 11) then begin
    printf,file_unit,"min_lat = -75"
    printf,file_unit,"max_lat = 75"
  endif else if (map eq 13) then begin
    printf,file_unit,"min_lat = -90 > (param1 - 90 + 15)"
    printf,file_unit,"max_lat =  90 < (param1 + 90 - 15)"
  endif else if (map eq 16) then begin
    printf,file_unit,"min_lat = -75"
    printf,file_unit,"max_lat = 90"
  endif else begin
    printf,file_unit,"min_lat = -90"
    printf,file_unit,"max_lat = 90"
  endelse
  printf,file_unit,"lat_spacing = 15"
endif else if (choice eq 2) then begin
  print,"Please enter the following quantities in degrees.'
  read,"  minimum longitude:",min_lon
  printf,file_unit,"min_lon = ",min_lon
  read,"  maximum longitude:",max_lon
  printf,file_unit,"max_lon = ",max_lon
  read,"  longitude spacing:",lon_spacing
  printf,file_unit,"lon_spacing = ",lon_spacing
  read,"  minimum latitude:",min_lat
  printf,file_unit,"min_lat = ",min_lat
  read,"  maximum latitude:",max_lat
  printf,file_unit,"max_lat = ",max_lat
  read,"  latitude spacing:",lat_spacing
  printf,file_unit,"lat_spacing = ",lat_spacing
endif

printf,file_unit,""
printf,file_unit,$
"; Based on the ranges for latitude and longitude, as well as their spacing,"
printf,file_unit,$
"; generate the latitude and longitude arrays."
printf,file_unit,"num_lon = long((max_lon - min_lon)/lon_spacing) + 1"
printf,file_unit,"lon = dindgen(num_lon)*lon_spacing + min_lon"
printf,file_unit,"num_lat = long((max_lat - min_lat)/lat_spacing) + 1"
printf,file_unit,"lat = dindgen(num_lat)*lat_spacing + min_lat"
printf,file_unit,"longitude = dblarr(num_lon,num_lat)"
printf,file_unit,"for i = 0,num_lat - 1 do longitude[*,i] = lon"
printf,file_unit,"latitude = dblarr(num_lon,num_lat)"
printf,file_unit,"for i = 0,num_lon - 1 do latitude[i,*] = lat"

printf,file_unit,""
printf,file_unit,";CONVERSION"

printf,file_unit,$
"; Convert the spherical coordinates into x-y coordinates by using wcssph2xy."
if (map lt 23) then begin
  if (n_elements(param1) eq 0) then begin
    if (n_elements(alpha) ne 0) then begin
      printf,file_unit,$
        "wcssph2xy,longitude,latitude,x,y,map,crval=[alpha,delta],$"
      printf,file_unit,"          longpole=longpole"
    endif else begin
      printf,file_unit,"wcssph2xy,longitude,latitude,x,y,map"
    endelse
  endif else if (n_elements(param2) eq 0) then begin
    if (n_elements(alpha) ne 0) then begin
      printf,file_unit,$
        "wcssph2xy,longitude,latitude,x,y,map,pv2=param1, $"
      printf,file_unit,"          crval=[alpha,delta],longpole=longpole"
    endif else begin
       printf,file_unit,"wcssph2xy,longitude,latitude,x,y,map,pv2=param1"
    endelse
  endif else begin
    if (n_elements(alpha) ne 0) then begin
      printf,file_unit,$
        "wcssph2xy,longitude,latitude,x,y,map,pv2=[param1,param2],$
      printf,file_unit,"          crval=[alpha,delta],longpole=longpole"
    endif else begin
      printf,file_unit,$
        "wcssph2xy,longitude,latitude,x,y,map,pv2=[param1,param2]"
    endelse
  endelse
endif else begin
  printf,file_unit,$
 "; The variable face must be declared with the same structure as latitude and"
  printf,file_unit,"; longitude before calling wcssph2xy." 
  printf,file_unit,"face = longitude - longitude"
  if (n_elements(alpha) ne 0) then begin
    printf,file_unit,$
      "wcssph2xy,longitude,latitude,x,y,map,face=face,crval=[alpha,delta], $
    printf,file_unit,"          longpole=longpole"
  endif else begin
    printf,file_unit,"wcssph2xy,longitude,latitude,x,y,map,face=face"
  endelse
endelse
printf,file_unit,""

printf,file_unit,";PLOTTING"
printf,file_unit,$
"; all maps have x increasing to the left, so switch this"
printf,file_unit,"xx = -x"
printf,file_unit,""

printf,file_unit,";LABELS"
printf,file_unit,$
"; The arrays lon_index and lat_index contain the indices for the latitude"
printf,file_unit,$
"; and longitude labels.  Labels occur every 30 degrees unless 30 doesn't"
printf,file_unit,$
"; divide into any of the latitude and longitude values evenly.  In this case,"
printf,file_unit,$
"; all latitude and longitude lines are labeled."
printf,file_unit,$
  "lon_index = where(long(longitude[*,0])/30 eq longitude[*,0]/30.)"
printf,file_unit,$
  "lat_index = where(long(latitude[0,*])/30 eq latitude[0,*]/30.)"
printf,file_unit,$
  "if (lat_index[0] eq -1) then lat_index = indgen(n_elements(latitude[0,*]))"
printf,file_unit,$
  "if (lon_index[0] eq -1) then lon_index = indgen(n_elements(longitude[*,0]))"

printf,file_unit,""

if (option lt 3) then begin
 if (n_elements(param2) eq 1) then wcssph2xy_plot,file_unit,map,param1,param2 $
  else if (n_elements(param1) eq 1) then wcssph2xy_plot,file_unit,map,param1 $
  else wcssph2xy_plot,file_unit,map

  if (option eq 2) then begin
    printf,file_unit,"key = ''"
    printf,file_unit,"read,'Press return to continue',key"

    if (n_elements(param2) eq 1) then $
       inversion_error,file_unit,map,param1,param2 $
    else if (n_elements(param1) eq 1) then $
       inversion_error,file_unit,map,param1 $
    else inversion_error,file_unit,map
  endif
endif else begin
  if (n_elements(param2) eq 1) then wcs_rot,file_unit,map,param1,param2 $
  else if (n_elements(param1) eq 1) then wcs_rot,file_unit,map,param1 $
  else wcs_rot,file_unit,map
endelse

printf,file_unit,"end"
close,file_unit
print,$
"The commands needed to execute what you are about to see can be executed"
print,"interactively, by typing ",strmid(file_name,0,strlen(file_name)-3)
print,""
command = strmid(file_name,0,strlen(file_name) - 4)
r = execute(command)
endrep until (option eq 4)
end