pro allsky_mag1,image,x,y,xc,yc,rad,maxzd,rot,h0,lat0,lon0,hl,hr,min,sec,xm,ym ; Procedure allsky_mag converts all-sky images at various ground stations ; from a horizon (az/el) system to a PACE geomagnetic coordinate system (Baker ; and Wing, JGR 94, 9139, 1989). ; ; Mapping is accomplished by producing geomagnetic look-up arrays that relate ; every point in the all-sky field of view with a corresponding magnetic ; latitude and magnetic local time. These look-up arrays will be used in ; conjunction with the program `polarplot3.pro' to display several allsky ; images on a common geomagnetic grid. Input variables that describe various ; image attributes are described as follows. ; ; Variable `image' is the allsky image array which is assumed to be a binary ; array. It is assumed that images have been photometrically processed, that ; is, the data number associated with all PIXELs can be ultimately related to ; Rayleighs. ; ; Variables `x' and `y' specify the dimensions of the all-sky image array, ; `xc' and `yc' are the PIXEL coordinates for the center field of view of the ; image, and `rad' is the PIXEL radius to the edge of the field of view. ; Image arrays are assumed to be in the normal `looking up' orientation, with ; east counterclockwise from north. Variable `rot' specifies the rotation of ; the CCD image plane y-axis away from geographic north and is given by the ; azimuth of the central column which passes through the center of the field ; of view (measured at the top of the image array). As an example, images ; oriented with geographic east at the top of the array would use a `rot' of ; 90 degrees. ; ; It is assumed that input image arrays have a 1:1 aspect ratio (the all-sky ; portion of such images arrays look like circles, rather than ovals) and that ; any effects from lens distortion have been taken out (the PIXEL radius from ; the image center is linearly related to true elevation angle). The zenith ; distance associated with edge of the field of view is given in variable ; `maxzd', and `maxzd' should be close to 90. degrees for truly all-sky images. ; ; Variables `h0', `lat0' and `lon0' are the altitude, geodetic latitude and ; geographic longitude for the allsky camera, respectively. Note that altitude ; should be given in kilometers and that west longitude should be specified ; as a negative number. ; ; Variable `h1' is supplied as the centroid height of the emission layer to be ; mapped. Typical values are 250 km for 6300 A and 110 km for 4278 A. ; Variables `hr', `min' and `sec' are the image UT times. ; ; xm and ym are output coordinate arrays that tell the plotting procedure ; (plot_image.pro) where to place each element of the image on a polar ; geomagnetic grid. This grid is oriented with the north geomagnetic pole at ; the center, 70 degrees magnetic latitude at the edge, and magnetic midnight ; at the bottom. ; ; ; R. Doe 8/22/94 SRI International ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Dimension arrays and define constants mlat = fltarr(x,y) mlt = fltarr(x,y) xm = fltarr(x,y) ym = fltarr(x,y) a = 6378.16 ab2 = 1.0067397 ep2 = 0.0067397 degrad= 180./!pi buf = fltarr(2) ; Load previously calculated geomagnetic MLAT and MLT lookup arrays ; These 730 by 798 floating point arrays relate geographic longitude (x-axis) ; and geographic latitude (y-axis) to PACE modeled MLAT and MLT. ; These lookup arrays have been calculated for 0 UT and 400 km. loadf,'mlat0.out',730,798,0,0,mlat0 loadf,'mlt0.out',730,798,0,0,mlt0 ; Transform the MLT array to correspond to the current UT time UT = float(hr) + float(min)/60. + float(sec)/3600. fac = !pi/12. sum = 0.0 param = [0.09, 0.14,-0.09, 0.04, 0.02, 0.02,-0.02, $ 0.00, 0.01, 0.01, 0.00, 0.00, 0.00, 0.00, $ 0.00, 0.00, 0.00, 0.00, 0.00] for n = 1,8 do begin sum = sum + param((2*n)-1)*sin(n*(UT*fac)) + $ param((2*n)) *cos(n*(UT*fac)) endfor correction = param(0) + sum mlt0 = mlt0 + UT + correction ind0 = where(mlt0 ge 24, count) if count gt 0 then mlt0(ind0) = mlt0(ind0) - 24.0 ; Calculate geocentric radius at transmitter and at emission layer cos_gdlat = cos(float(lat0)/degrad ) sin_gdlat = sin(float(lat0)/degrad ) sl2 = (sin_gdlat)^2. cl2 = (ab2 * cos_gdlat)^2. sinbet = (sin_gdlat)/sqrt(sl2 + cl2) buf(0:1) = [(sinbet)^2.,1.] sb2 = min(buf) cosbet = sqrt(1. - sb2) rgeoid = a/sqrt(1. + ep2 * sb2) xx = rgeoid * cosbet + float(h0) * cos_gdlat yy = rgeoid * sinbet + float(h0) * sin_gdlat t_rad = sqrt(xx^2. + yy^2.) l_rad = t_rad + float(hl) ; Determine the elevation angle, azimuth angle and range associated with each ; point in the image array, convert to geodetic latitude, longitude and ; altitude (lat1, lon1, alt1), and then convert to magnetic latitude and ; magnetic local time (mlat0 and mlt0) using the previously defined look-up ; arrays. for i = 0,x-1 do begin print,i,x-1 for j = 0,y-1 do begin xxx = float(i) - float(xc) yyy = float(j) - float(yc) r = sqrt(xxx^2. + yyy^2.) if r gt float(rad) then goto, continue0 theta = atan(yyy,xxx+.00001)*(180./!pi) if theta lt 0. then theta = 360. + theta azimuth = theta - 90. + float(rot) if azimuth lt 0. then azimuth = 360. + azimuth elevation = 90. - r * (float(maxzd)/float(rad)) if elevation lt 0. then elevation = 0. az0 = azimuth el0 = elevation beta = 180.-(90.+el0)-asin((t_rad/l_rad)*sin((el0+90.)/degrad))*degrad range = sqrt(t_rad^2.+l_rad^2.-2.*t_rad*l_rad*cos(beta/degrad)) elaz1,float(lat0),float(lon0),float(h0),el0,az0,range,lat1,lon1,alt1 if lon1 lt 0. then lon1 = 360. + lon1 ; Transform lon1 and lat1 to corresponding look-up array addresses xxxx = lon1 * (729./360.) yyyy = (90. - lat1) * (797./47.5) if (xxxx le 729.) and (yyyy le 797.) then begin mlat1 = mlat0(xxxx,yyyy) mlt1 = mlt0(xxxx,yyyy) endif mlat(i,j) = mlat1 mlt(i,j) = mlt1 continue0: endfor endfor radius = (90. - mlat)*(799./40.) theta = (mlt * (360./24.)) - 90. xm = 399. + radius * cos(theta * (!pi/180.)) ym = 399. + radius * sin(theta * (!pi/180.)) indx = where((xm lt 0.) or (xm gt 799.),count) if count gt 0 then xm(indx) = 0. indy = where((ym lt 0.) or (ym gt 799.),count) if count gt 0 then ym(indy) = 0. end