;+ ; NAME: ; swand ; ; PURPOSE: ; ; This file contains all the SWAN software required to read and analyse ; level 0 SWAN fits files within IDL. ; ; This widget program illustrates the use of this software. It shows a ; complete analysis chain for full sky observations. ; ; To run the program just type "swand" at the IDL prompt. Then select a ; fullsky from the menu and then proceed. ; ; You need to have the IDL astronomy user's library installed ; on your system. ; ; A number of "IDL save" files are required (see the SWAN user's ; resource page): ; - the orbit of SoHO is defined in orbitsohoidl.dat, ; - rollangle.dat define the nominal roll angle of SoHO, ; - MasqueBaf2.dat contains two stars masks obtained from ; the SU+Z pixel 26 counting rate (upper thresholds 10 ; and 12 cps). This pixel has a BaF2 window which cuts off ; Lyman alpha emission. It is used to track stellar light. ; - offset-Zidl.dat, offset+Zidl.dat are used to correct ; misalignments of the optical path in both sensors. ; ; The file orbitsohoidl.dat is updated regularly. ; ; Since the summer of 2003, SOHO has to roll around the Ox axis ; (sun pointing axis) during half of its halo orbit in order to ; maximize the data rate of the high gain antenna. The nominal ; roll angle switches between 0 and 180 degrees. The values are ; defined in the file rollangle.dat. Its content has to be ; updated every three months. ; ; CATEGORY: ; Widget ; ; CALLING SEQUENCE: ; swand ; ; INPUTS: ; The paths of FITS and "IDL save" files . ; ; KEYWORD PARAMETERS: ; None ; ; OUTPUTS: ; The user can save a copy of the full sky image (gif) and/or ; save the following data within a IDL save file (bin): ; - tfskz, tfsk_z: begining of the SU+Z and SU-Z observations ; (in number of days since 01/01/1995), ; - cmean: fullsky image, cmean(*,*,0) contains the data, ; cmean(*,*,1) contains the relative statistical error ; - tnord, tsud: average time of observation for each ligne ; of sight from SU+Z (tnord) and SU-Z (tsud) ; - titreplot: some title for plotting, ; - flatsuz, flatsu_z: flatfield corrections values, errors ; and reduced chisquare, ; - coefns: intercalibration factor between -Z (south) and +Z ; (north) sensors (values, errors and reduced chisquare), ; - cal: absolute calibration factor, if determined, with ; its error. It applies to the count rate of the +Z (north) ; sensor. ; ; RESTRICTIONS: ; The program is intended to work only on full sky observations. ; ; The analysis of other modes (eg comet observations) requires ; a prior determination of the flatfield corrections from full ; sky observations taken during the same period and in the same ; high voltage conditions. ; ; PROCEDURES/FUNCTIONS: ; - read_from_fits: read data from SWAN FITS files, ; - invdat: given the number of days since 1995/01/01, extract ; the year, the month and day in the month, ; - jourdec2txt: convert a number of days since 1995/01/01 into ; a readable string, ; - cleanid: apply selections on telemetry packets, ; - avefixmot: apply corrections on motors positions, ; - checkoutinn: correct any motors inversion for SU-Z, ; - trouve3d: linear interpolations on three arrays, ; - motoffset: apply corrections on motors positions, ; - reflection: computes the direction of a light ray reflected ; on a mirror, ; - su2soho27: computes the transformation from motor angles and ; pixel position into SOHO reference frame, ; - soho2ecli: computes the transformation from SOHO coordinates ; to ecliptic coordinates (J2000), ; - sohoroll2ecli: same as previous. Must be used in case of non ; zero roll angle, ; - convert3to2deg, convert3to2deg0: conversion from cartesian to ; spherical coordinates (in degrees), ; - ang2vec: angle between two vectors, ; - fpr: robust fit used for flatfield and N/S corrections ; determinations, ; - keeparea: forbidden area definition (avoid light scattered by ; the S/C), ; - flficomp: flatfield fits, ; - ecli2equa: ecliptic to Earth Equatorial coordinates (J2000) ; (or inverse transformation), ; - ecli2wind: ecliptic to wind-Oy coordinates (or inverse ; transformation), ; - equa2gala: earth Equatorial coordinates to Galactic ; coordinates (or inverse transformation), ; - mkim: image calculation, ; - htd2a: conversion between analogique and digital high voltage ; values, ; - plotimage: 2D representation of an image, ; - analysis: the main analysis routine, ; - SwandEventHandler: event handler, ; - swand: widget setup. ; ; MODIFICATION HISTORY: ; Written by: Cyril Pennanech, Eric Quemerais, Stephane Ferron ; Test on IDL Version 6.1 (linux x86 m32, windows XP). ; ; 2007-12-06 SF idlsave file name and paths mods. ;$Id: swand.pro,v 1.9 2005/09/26 14:49:54 iws Exp $ ; ;- pro read_from_fits, file, widget_io = widget_io common vardet, packet, secondes, detsuz, detsu_z, insuz, outsuz, insu_z, $ outsu_z, suzhvps, su_zhvps, suzhcell, su_zhcell, tiz, ti_z, $ outmotz, innmotz, outmot_z, innmot_z, nbre, id, iprim, $ sec_det common header, hearder, ha, hb, hc, hd, he, hf, hg, hh, hi, hj, hk, hl, $ hm, hn, ho, hp, hq, hr, hs, ht if keyword_set(widget_io) then widget_control, widget_io, set_value = 'Opening ' + file + '.', /append else print, 'read_from_fits: opening ' + file + '.' ;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Read header in fits file ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;; header = headfits(file) packet = fxpar(header,'SW_NPACK') senseur = fxpar(header,'DETECTOR') case 1 of strpos(senseur,'BOTH') gt -1 : z = 2 strpos(senseur,'SU+Z') gt -1 : z = 1 strpos(senseur,'SU-Z') gt -1 : z = 0 strpos(senseur,'No DATA') gt -1 : z = -1 endcase data_no = 1 if z eq -1 then begin if keyword_set(widget_io) then widget_control, widget_io, set_value = 'No data from this fits file.', /append else print, 'No data from this fits file.' goto, theend endif ;;;;;;;;;;;;;;;;;; ; Read SU+Z data ; ;;;;;;;;;;;;;;;;;; if z eq 1 or z eq 2 then begin detsuz = readfits(file,ha,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 outsuz = readfits(file,hb,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 insuz = readfits(file,hc,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 outmotz = readfits(file,hd,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 innmotz = readfits(file,he,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 suzhcell = readfits(file,hf,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 suzhvps = readfits(file,hg,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 tiz = readfits(file,hh,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 endif else begin detsuz = -1 outsuz = -1 insuz = -1 outmotz = -1 innmotz = -1 suzhcell = -1 suzhvps = -1 tiz = -1 endelse ;;;;;;;;;;;;;;;;;; ; Read SU-Z data ; ;;;;;;;;;;;;;;;;;; if z eq -1 or z eq 2 then begin detsu_z = readfits(file,hi,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 outsu_z = readfits(file,hj,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 insu_z = readfits(file,hk,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 outmot_z = readfits(file,hl,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 innmot_z = readfits(file,hm,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 su_zhcell = readfits(file,hn,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 su_zhvps = readfits(file,ho,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 ti_z = readfits(file,hp,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 endif else begin detsu_z = -1 outsu_z = -1 insu_z = -1 outmot_z = -1 innmot_z = -1 su_zhcell = -1 su_zhvps = -1 ti_z = -1 endelse ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Read Common data For both Sensor ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; secondes = readfits(file,hq,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 id = readfits(file,hr,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 iprim = readfits(file,hs,EXTEN_NO = data_no,/SILENT) data_no = data_no + 1 sec_det = readfits(file,ht,EXTEN_NO = data_no,/SILENT) theend: end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function invdat, date yr = 0 mo = 0 ld = 2 fr = date - floor(date) dn = date - fr - 365. while dn ge 365. * (yr + 1) + floor(yr / 4.) + 1 do yr = yr + 1 dn = dn - 365. * yr - floor((yr - 1) / 4.) - 1 if (yr mod 4) eq 0 then ld = 1 while dn ge 30. * (mo + 1) + floor(mo / 2.) + 1 - ld * (mo ge 1) - 30 * (mo ge 7) do mo = mo + 1 dn = dn - 30. * mo - floor((mo - 1) / 2.) - 1 + ld * (mo gt 1) + 30 * (mo gt 7) if mo lt 8 then mo = mo + 1 return, [1996. + yr, mo, dn + 1, 24 * fr] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function jourdec2txt, jourdec fdate = invdat(jourdec) return, string(fdate(0),format = '(i4)') + '/' + string(fdate(1), format = '(i02)') + '/' + string(fdate(2), format = '(i02)') end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro cleanid, su, id, jourdec, jourdec1, outsuz, outmotz, insuz, innmotz, detsuz, suzhcell, suzhvps, tiz, widget_io = widget_io nbre = n_elements(outsuz) ht1 = where(suzhvps lt 1600., nht1) if nht1 ge 1 then begin for i = 0, nht1 - 1 do begin jj = ht1(i) if keyword_set(widget_io) then $ widget_control, widget_io, set_value = 'SU' + (su eq 1 ? '+' : '-') + 'Z HV OFF at Packet ' + string(jj, format = '(i5)') + '.', /append else $ print, 'SU' + (su eq 1 ? '+' : '-') + 'Z HV OFF at Packet ' + string(jj, format = '(i5)') + '.' suzhvps(jj) = 0. if (jj - 1 ge 0) then suzhvps(jj - 1) = 0. if (jj + 1 lt nbre) then suzhvps(jj + 1) = 0. endfor endif id1 = where(suzhvps gt 1600. and id eq 1, nid1) if nid1 gt 1 then begin nid1 = nid1 - 1 ; skip the first ... id1 = id1(1:nid1) ; and the last packet. outsuz = outsuz(id1) insuz = insuz(id1) outmotz = outmotz(id1) innmotz = innmotz(id1) dummy = fltarr(nid1,27) for i = 0, nid1 - 1 do for j = 0, 26 do dummy(i,j) = detsuz(id1(i),j) detsuz = dummy jourdec1 = jourdec(id1) suzhcell = suzhcell(id1) suzhvps = suzhvps(id1) tiz = tiz(id1) nbre1 = n_elements(id1) endif else begin nbre1 = 0 outsuz = 0. insuz = 0. endelse end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro avefixmot, su, titrav, outsuz, outmotz, plot = plot outsuzold = outsuz nbre1 = n_elements(outsuz) n0 = 0 newbin0: n1 = n0 newbin: n1 = n1 + 1 if (outmotz(n1) eq outmotz(n0) and n1 lt nbre1 - 1) then goto, newbin n1 = n1 - 1 bb = where(abs(outsuz(n0:n1) - outmotz(n0:n1)) lt 160., nbe) if nbe gt 1 then begin bb = n0 + bb amoy = total(float(outsuz(bb))) / float(nbe) diffmax = max(abs(outsuz(bb) - amoy)) if diffmax lt 10 then outsuz(bb) = amoy endif n0 = n1 + 1 if n0 lt nbre1 - 1 then goto, newbin0 if keyword_set(plot) then plot, outsuz - outsuzold, xtitle = 'Packet number', ytitle = 'Motor correction', title = titrav + (su eq -1 ? ' SU-Z' : ' SU+Z') end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro checkoutinn, outsu_z, outmot_z, insu_z, innmot_z, widget_io = widget_io eval1 = total(abs(outsu_z - outmot_z)) / n_elements(outsu_z) eval2 = total(abs(outsu_z - innmot_z)) / n_elements(outsu_z) if eval1 gt eval2 then begin dum = outmot_z outmot_z = innmot_z innmot_z = dum if keyword_set(widget_io) then widget_control, widget_io, set_value = ' The innmot_z and outmot_z variables have been inverted.', /append else $ print,'checkoutinn: the innmot_z and outmot_z variables have been inverted.' endif end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function trouve3d, date, temps, xcystem, ycystem, zcystem ntem = n_elements(temps) - 1 if date le temps(0) then return, [xcystem(0), ycystem(0), zcystem(0)] if date ge temps(ntem) then return, [xcystem(ntem), ycystem(ntem), zcystem(ntem)] bon = where(temps le date, nbon) n1 = max(bon) delt = (date - temps(n1)) / (temps(n1 + 1) - temps(n1)) x0 = xcystem(n1) * (1. - delt) + xcystem(n1 + 1) * delt y0 = ycystem(n1) * (1. - delt) + ycystem(n1 + 1) * delt z0 = zcystem(n1) * (1. - delt) + zcystem(n1 + 1) * delt return, [x0, y0, z0] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function motoffset, suz, outin, innin common offset, outoffz, innoffz, outoff_z, innoff_z, xout, xinn, nout, ninn outout = outin innout = innin outmin = -90. & outdel = 10. & outmax = 90. innmin = -90. & inndel = 10. & innmax = 90. if suz eq 1 then begin outoff = outoffz innoff = innoffz endif else begin outoff = outoff_z innoff = innoff_z endelse out = outout & inn = innout if abs(out) ge 5. then begin jout = (out - outmin) / outdel outi = fix(jout) & outs = fix(jout) + 1 if outi lt 0 then outi = 0 if outi gt nout - 1 then outi = nout - 1 if outs lt 0 then outs = 0 if outs gt nout - 1 then outs = nout - 1 dout = jout - 1. * outi endif else begin if out lt 0. then begin outi = 8 & outs = 8 & dout = 0. endif else begin outi = 9 & outs = 9 & dout = 0. endelse endelse jinn = (inn - innmin) / inndel inni = fix(jinn) & inns = fix(jinn) + 1 if inni lt 0 then inni = 0 if inni gt ninn - 1 then inni = ninn - 1 if inns lt 0 then inns = 0 if inns gt ninn - 1 then inns = ninn - 1 dinn = jinn - 1. * inni out1 = outoff(outi,inni) * (1. - dout) * (1. - dinn) + $ outoff(outs,inni) * dout * (1. - dinn) + $ outoff(outi,inns) * (1. - dout) * dinn + $ outoff(outs,inns) * dout * dinn inn1 = innoff(outi,inni) * (1. - dout) * (1. - dinn) + $ innoff(outs,inni) * dout * (1. - dinn) + $ innoff(outi,inns) * (1. - dout) * dinn + $ innoff(outs,inns) * dout * dinn outout = out + out1 / 40. innout = inn + inn1 / 40. return, [outout, innout] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function reflection, normale, vector_in uu = normale(0) vv = normale(1) ww = normale(2) out1 = (2. * uu * uu - 1.) * vector_in(0) + 2. * uu * vv * vector_in(1) + 2. * uu * ww * vector_in(2) out2 = 2. * uu * vv * vector_in(0) + (2. * vv * vv - 1.) * vector_in(1) + 2. * vv * ww * vector_in(2) out3 = 2. * uu * ww * vector_in(0) + 2. * ww * vv * vector_in(1) + (2. * ww * ww - 1.) * vector_in(2) return, [out1, out2, out3] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function su2soho27, outrot, inrot, suz sortie = fltarr(3,27) focale = 108. ; millimeters theta = inrot * !pi / 180. aphi = outrot * !pi / 180. normal_1=[cos(theta), sin(theta), -1.] / sqrt(2.) normal_2=[-cos(theta) + sin(theta) * sin(aphi), $ -sin(theta) - cos(theta) * sin(aphi), $ cos(aphi)] / sqrt(2.) xpixel = 1.78 * float(2 - indgen(25) / 5) ypixel = 1.78 * ((findgen(25) mod 5.) - 2.) longueur=(1. / sqrt(xpixel^2 + ypixel^2 + focale^2)) for jj = 0, 24 do begin los0 = [xpixel(jj), ypixel(jj), focale] * longueur(jj) los1 = reflection(normal_1,-los0) los2 = reflection(normal_2,-los1) norme = sqrt(total(los2 * los2)) sortie(0,jj) = los2(0) / norme sortie(1,jj) = suz * los2(1) / norme sortie(2,jj) = suz * los2(2) / norme endfor jj = 25 xpix = 0. if suz eq 1 then ypix = -(3. * 1.78 + 1.8) else ypix = 3. * 1.78 + 1.8 longur = (1. / sqrt(xpix^2 + ypix^2 + focale^2)) los0 = [xpix, ypix, focale] * longur los1 = reflection(normal_1,-los0) los2 = reflection(normal_2,-los1) norme = sqrt(total(los2 * los2)) sortie(0,jj) = los2(0) / norme sortie(1,jj) = suz * los2(1) / norme sortie(2,jj) = suz * los2(2) / norme jj = 26 if suz eq 1 then xpix = 3. * 1.78 + 0.8 else xpix = 0. if suz eq 1 then ypix = 0. else ypix = -3. * 1.78 + 1.8 longur = (1. / sqrt(xpix^2 + ypix^2 + focale^2)) los0 = [xpix, ypix, focale] * longur los1 = reflection(normal_1,-los0) los2 = reflection(normal_2,-los1) norme = sqrt(total(los2 * los2)) sortie(0,jj) = los2(0) / norme sortie(1,jj) = suz * los2(1) / norme sortie(2,jj) = suz * los2(2) / norme return, sortie end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function soho2ecli, sc_posit, coord_in, inv dso = sqrt(total(sc_posit^2)) i1s = [-sc_posit(0), $ -sc_posit(1), $ -sc_posit(2)] / dso k1s = (-i1s(2) * i1s + [0., 0., 1.]) / sqrt(1. - i1s(2)^2) j1s = [k1s(1) * i1s(2) - k1s(2) * i1s(1), $ k1s(2) * i1s(0) - k1s(0) * i1s(2), $ k1s(0) * i1s(1) - k1s(1) * i1s(0)] omega = 1.31982 incli = 0.126536 ksun = [sin(omega) * sin(incli), -cos(omega) * sin(incli), cos(incli)] beta = ksun(0) * j1s(0) + ksun(1) * j1s(1) + ksun(2) * j1s(2) gamma = ksun(0) * k1s(0) + ksun(1) * k1s(1) + ksun(2) * k1s(2) ks = beta * j1s + gamma * k1s norme = sqrt(total(ks * ks)) ks = ks / norme is = i1s js = [ks(1) * is(2) - ks(2) * is(1), $ ks(2) * is(0) - ks(0) * is(2), $ ks(0) * is(1) - ks(1) * is(0)] if inv eq 1 then begin x_out = coord_in(0,*) * is(0) + coord_in(1,*) * js(0) + coord_in(2,*) * ks(0) y_out = coord_in(0,*) * is(1) + coord_in(1,*) * js(1) + coord_in(2,*) * ks(1) z_out = coord_in(0,*) * is(2) + coord_in(1,*) * js(2) + coord_in(2,*) * ks(2) endif else begin x_out = coord_in(0,*) * is(0) + coord_in(1,*) * is(1) + coord_in(2,*) * is(2) y_out = coord_in(0,*) * js(0) + coord_in(1,*) * js(1) + coord_in(2,*) * js(2) z_out = coord_in(0,*) * ks(0) + coord_in(1,*) * ks(1) + coord_in(2,*) * ks(2) endelse return, [x_out, y_out, z_out] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function sohoroll2ecli, sc_posit, coord_in, inv, angroll dso = sqrt(total(sc_posit^2)) i1s = [-sc_posit(0), $ -sc_posit(1), $ -sc_posit(2)] / dso k1s = (-i1s(2) * i1s + [0., 0., 1.]) / sqrt(1. - i1s(2)^2) j1s = [k1s(1) * i1s(2) - k1s(2) * i1s(1), $ k1s(2) * i1s(0) - k1s(0) * i1s(2), $ k1s(0) * i1s(1) - k1s(1) * i1s(0)] omega = 1.31982 incli = 0.126536 ksun = [sin(omega) * sin(incli), -cos(omega) * sin(incli), cos(incli)] beta = ksun(0) * j1s(0) + ksun(1) * j1s(1) + ksun(2) * j1s(2) gamma = ksun(0) * k1s(0) + ksun(1) * k1s(1) + ksun(2) * k1s(2) ks = beta * j1s + gamma * k1s norme = sqrt(total(ks*ks)) ks = ks / norme is = i1s js = [ks(1) * is(2) - ks(2) * is(1), $ ks(2) * is(0) - ks(0) * is(2), $ ks(0) * is(1) - ks(1) * is(0)] irs = is jrs = cos(angroll * !dtor) * js + sin(angroll * !dtor) * ks krs = -sin(angroll * !dtor) * js + cos(angroll * !dtor) * ks if inv eq 1 then begin x_out = coord_in(0,*) * irs(0) + coord_in(1,*) * jrs(0) + coord_in(2,*) * krs(0) y_out = coord_in(0,*) * irs(1) + coord_in(1,*) * jrs(1) + coord_in(2,*) * krs(1) z_out = coord_in(0,*) * irs(2) + coord_in(1,*) * jrs(2) + coord_in(2,*) * krs(2) endif else begin x_out = coord_in(0,*) * irs(0) + coord_in(1,*) * irs(1) + coord_in(2,*) * irs(2) y_out = coord_in(0,*) * jrs(0) + coord_in(1,*) * jrs(1) + coord_in(2,*) * jrs(2) z_out = coord_in(0,*) * krs(0) + coord_in(1,*) * krs(1) + coord_in(2,*) * krs(2) endelse return, [x_out, y_out, z_out] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function convert3to2deg, coordinate norme = sqrt(coordinate(0,*)^2 + coordinate(1,*)^2 + coordinate(2,*)^2) xecli = coordinate(0,*) / norme yecli = coordinate(1,*) / norme zecli = coordinate(2,*) / norme lon = atan(yecli,xecli) * 180. / !pi corr = where(lon lt 0., ncor) if ncor gt 0 then lon(corr) = lon(corr) + 360. lat = asin(zecli) * 180. / !pi return, [lon, lat] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function convert3to2deg0, coordinate norme = sqrt(coordinate(0,*)^2 + coordinate(1,*)^2 + coordinate(2,*)^2) xecli = coordinate(0,*) / norme yecli = coordinate(1,*) / norme zecli = coordinate(2,*) / norme lon = atan(yecli,xecli) * 180. / !pi lat = asin(zecli) * 180. / !pi return, [lon, lat] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function ang2vec, vector1, vector2 t1 = size(vector1) t2 = size(vector2) ang = 0. if (t1(1) eq 3 and t2(1) eq 3) then begin n1 = sqrt(total(vector1 * vector1)) n2 = sqrt(total(vector2 * vector2)) pdsc = total(vector1 * vector2) / n1 / n2 if abs(pdsc) lt 1. then ang = acos(pdsc) / !dtor else begin if pdsc ge 1. then ang = 0. if pdsc le -1. then ang = 180. endelse endif return, ang end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function fpr, vx, ex, vy, ey, ns kp = where(vx ne 0., np) if np le 1 then begin print, 'fpr: not enough non zero values in vx!' return, -1 endif vr = vy(kp) / vx(kp) er = vr * sqrt(ex(kp)^2 + ey(kp)^2) df = float(np - 1) s0 = total(1. / er^2) s1 = total(vr / er^2) s2 = total(vr^2 / er^2) va = s1 / s0 ea = sqrt(df / s0) cn = (s2 - s1^2 / s0) / df ; print, 'fpr: first pass, (chi2/n.d.f. = ' + string(cn, format = '(f)') + '), r = ' + string(va, format = '(f7.3)') + ' +/- ' + string(ea, format = '(f)') + '.' kp = where(abs(vr - va) le ns * ea, np) if np le 1 then begin print, 'fpr: not enough non zero values in vr!' return, -1 endif df = float(np - 1) s0 = total(1. / er(kp)^2) s1 = total(vr / er(kp)^2) s2 = total(vr^2 / er(kp)^2) va = s1 / s0 ea = sqrt(df / s0) cn = (s2 - s1^2 / s0) / df ; print, 'fpr: second pass, (chi2/n.d.f. = ' + string(cn, format = '(f)') + '), r = ' + string(va, format = '(f7.3)') + ' +/- ' + string(ea, format = '(f)') + '.' return, [va, ea, cn] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function keeparea, su, inn, out if su eq 1 then begin if (out lt 3000 and inn gt 7000) then return, 0 if (out ge 8200 and inn le 3200) then return, 0 if (out lt 3200 and inn gt 8400) then return, 0 if (out ge 8800 and inn le 5800) then return, 0 if (out ge 8500 and inn le 4800) then return, 0 if (out gt 8600 and inn gt 8400) then return, 0 if (out lt 3400 and inn lt 3500) then return, 0 if (out le 3000 and inn le 3800) then return, 0 if (out ge 9000 and inn ge 8200) then return, 0 if (out le 2800 and inn ge 6600) then return, 0 if (out le 3600 and inn ge 9400) then return, 0 endif else begin if (out lt 2800 and inn lt 3200) then return, 0 if (out gt 9200 and inn gt 8800) then return, 0 if (out gt 9500 and inn lt 5400) then return, 0 if (out lt 2500 and inn gt 6600) then return, 0 if (out lt 3000 and inn gt 7600) then return, 0 if (out gt 9000 and inn lt 3800) then return, 0 endelse return, 1 end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function flficomp, suz, countf, d2countf, innef, outef, tempf, masque common position, tempsoho, xx_soho, yy_soho, zz_soho, nombre common poterre, xterre_soho, yterre_soho, zterre_soho common rollangle, tfroll, rollval, optroll nlon = 360 lonmin = 0. londel = 1. nlat = 180 latmin = -90. latdel = 1. nlos = long(nlon) * long(nlat) pixbuf = fltarr(25,nlos,2) nsigm = 2.5 cwin = fltarr(25,2) flfiv = fltarr(25,3) coupoff = fltarr(25,nlon + 1,nlat,2) coupoff(0:24,nlon,0,1) = lonmin coupoff(0:24,nlon,1,1) = londel coupoff(0:24,nlon,2,1) = latmin coupoff(0:24,nlon,3,1) = latdel coupaux = fltarr(27,nlon + 1,nlat) ter_pos = trouve3d(tempf(0),tempsoho,xterre_soho,yterre_soho,zterre_soho) if optroll eq 1 then begin nrol = n_elements(rollval) if nrol lt 1 then optroll = 0 else begin nrl = where(tfroll le tempf(0),nnl) if nnl gt 0 then nrl = max(nrl) else nrl = 0 endelse endif for ll = 0, n_elements(tempf) - 1 do begin scposit = trouve3d(tempf(ll),tempsoho,xx_soho,yy_soho,zz_soho) out = outef(ll) inn = innef(ll) vout = motoffset(suz,out,inn) out = vout(0) inn = vout(1) tsout = fix(6000 + 40 * out) tsinn = fix(6000 + 40 * inn) if not keeparea(suz,tsinn,tsout) then goto, suivant visee1 = su2soho27(out,inn,suz) anrol = 0. if optroll eq 1 then begin suitrol: if nrl lt nrol - 1 then begin if tempf(ll) ge tfroll(nrl + 1) then begin nrl = nrl + 1 goto, suitrol endif endif anrol = rollval(nrl) visee3 = sohoroll2ecli(scposit,visee1,1,anrol) endif else visee3 = soho2ecli(scposit,visee1,1) visee1 = visee3 visee2 = convert3to2deg(visee1) for rr = 0, 24 do begin vtest = [visee1(0,rr), visee1(1,rr), visee1(2,rr)] vtes1 = ang2vec(vtest,ter_pos) if vtes1 lt 20. then goto, suivantrr jlon = (visee2(0,rr) - lonmin) / londel jlat = (visee2(1,rr) - latmin) / latdel if (jlon ge 0. and jlat gt 0. and jlon le nlon - 1 and jlat le nlat - 1) then begin lon = fix(jlon) lat = fix(jlat) coupoff(rr,lon,lat,0) = coupoff(rr,lon,lat,0) + countf(ll,rr) * masque(lon,lat) coupoff(rr,lon,lat,1) = coupoff(rr,lon,lat,1) + d2countf(ll,rr) * masque(lon,lat) coupaux(rr,lon,lat) = coupaux(rr,lon,lat) + 1. * masque(lon,lat) endif suivantrr: endfor suivant: endfor for rr = 0, 24 do begin for ii = 0, nlat - 1 do begin for jj = 0, nlon - 1 do begin if coupoff(rr,jj,ii,0) gt 0. then begin coupoff(rr,jj,ii,1) = sqrt(coupoff(rr,jj,ii,1)) / coupoff(rr,jj,ii,0) coupoff(rr,jj,ii,0) = coupoff(rr,jj,ii,0) / coupaux(rr,jj,ii) endif else coupoff(rr,jj,ii,*) = 0. if abs(latmin + (float(ii) + 0.5)) le 40. then coupoff(rr,jj,ii,0:1) = 0. pixbuf(0:24,long(ii) * long(nlon) + long(jj),0:1) = coupoff(0:24,jj,ii,0:1) endfor endfor tmpbuf = pixbuf(rr,*,0) histo = fltarr(1000) hmean = 0. hsigm = 0. for k = 0L, nlos - 1L do begin if ((tmpbuf(k) gt 0.5) and (tmpbuf(k) lt 1000.5)) then begin ix = fix(tmpbuf(k)) if (ix lt 1000) then histo(ix) = histo(ix) + 1. endif endfor ix = 0.5 + findgen(1000) hmean = total(ix * histo) / total(histo) hsigm = total(ix * ix * histo) / total(histo) hsigm = sqrt(hsigm - hmean ^ 2) m1 = min([max([fix(hmean - nsigm * hsigm), 0]),999]) m2 = min([max([fix(hmean + nsigm * hsigm), 0]),999]) if m1 gt 0 then histo(0:m1) = 0. if m2 lt 1000 then histo(m2:999) = 0. hmean = total(ix * histo) / total(histo) hsigm = total(ix * ix * histo) / total(histo) hsigm = sqrt(hsigm - hmean ^ 2) clean = where(tmpbuf gt hmean + nsigm * hsigm, nclean) if nclean gt 0 then pixbuf(rr,clean,0:1) = 0. clean = where(tmpbuf lt hmean - nsigm * hsigm, nclean) if nclean gt 0 then pixbuf(rr,clean,0:1) = 0. cwin(rr,0) = max([hmean - nsigm * hsigm, 0]) cwin(rr,1) = min([hmean + nsigm * hsigm, 1000]) endfor ref_val = pixbuf(12,*,0) ref_err = pixbuf(12,*,1) !p.multi = [0, 5, 5] for rr = 0, 24 do begin if rr eq 12 then begin flfiv(rr,0:2) = 1. !p.multi(0) = !p.multi(0) - 1 endif else begin cur_val = pixbuf(rr,*,0) cur_err = pixbuf(rr,*,1) comm = where(ref_err gt 0. and cur_err gt 0., npoint) if npoint le 1 then begin print, 'flficomp: no data for pixel ' + string(rr, format = '(i2)') + '.' flfiv(rr,0:2) = -1. !p.multi(0) = !p.multi(0) - 1 endif else begin flfi = fpr(cur_val(comm),cur_err(comm),ref_val(comm),ref_err(comm),2.) plot, ref_val(comm), cur_val(comm), psy = 3, xrange = [cwin(12,0), cwin(12,1)], yrange = [cwin(rr,0), cwin(rr,1)], xtitle = 'pixel 12', ytitle = 'pixel ' + string(rr, format = '(i2)') if n_elements(flfi) eq 3 then begin flfiv(rr,0) = flfi(0) flfiv(rr,1) = flfi(1) flfiv(rr,2) = flfi(2) endif else begin flfiv(rr,0) = 1. flfiv(rr,1) = 0. flfiv(rr,2) = 0. endelse pts = [min([cwin(12,0), cwin(rr,0)]), max([cwin(12,1), cwin(rr,1)])] oplot, pts, pts / flfi(0) oplot, pts, pts, linestyle = 2 endelse endelse endfor return, flfiv end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function ecli2equa, coord_in, inv is = [1., 0., 0.] js = [0., 0.917408, 0.397949] ks = [0., -0.397949, 0.917408] if inv eq 1 then begin x_out = coord_in(0,*) * is(0) + coord_in(1,*) * js(0) + coord_in(2,*) * ks(0) y_out = coord_in(0,*) * is(1) + coord_in(1,*) * js(1) + coord_in(2,*) * ks(1) z_out = coord_in(0,*) * is(2) + coord_in(1,*) * js(2) + coord_in(2,*) * ks(2) endif else begin x_out = coord_in(0,*) * is(0) + coord_in(1,*) * is(1) + coord_in(2,*) * is(2) y_out = coord_in(0,*) * js(0) + coord_in(1,*) * js(1) + coord_in(2,*) * js(2) z_out = coord_in(0,*) * ks(0) + coord_in(1,*) * ks(1) + coord_in(2,*) * ks(2) endelse return, [x_out, y_out, z_out] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function ecli2wind, coord_in, inv is = [ -0.961262, -0.273279, 0.0359779] js = [ 0.275637, -0.953038, 0.125470] ks = [ 0., 0.130526, 0.991445] if inv eq 1 then begin x_out = coord_in(0,*) * is(1) + coord_in(1,*) * js(1) + coord_in(2,*) * ks(1) y_out = -coord_in(0,*) * is(0) - coord_in(1,*) * js(0) - coord_in(2,*) * ks(0) z_out = coord_in(0,*) * is(2) + coord_in(1,*) * js(2) + coord_in(2,*) * ks(2) endif else begin x_out = coord_in(0,*) * is(1) - coord_in(1,*) * is(0) + coord_in(2,*) * is(2) y_out = coord_in(0,*) * js(1) - coord_in(1,*) * js(0) + coord_in(2,*) * js(2) z_out = coord_in(0,*) * ks(1) - coord_in(1,*) * ks(0) + coord_in(2,*) * ks(2) endelse return, [x_out, y_out, z_out] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function equa2gala, coord_in, inv is = [-0.0671515, -0.872716, -0.483588] js = [0.493086, -0.450385, 0.744325] ks = [-0.867436, -0.189132, 0.460200] if inv eq -1 then begin x_out = coord_in(0,*) * is(0) + coord_in(1,*) * js(0) + coord_in(2,*) * ks(0) y_out = coord_in(0,*) * is(1) + coord_in(1,*) * js(1) + coord_in(2,*) * ks(1) z_out = coord_in(0,*) * is(2) + coord_in(1,*) * js(2) + coord_in(2,*) * ks(2) endif else begin x_out = coord_in(0,*) * is(0) + coord_in(1,*) * is(1) + coord_in(2,*) * is(2) y_out = coord_in(0,*) * js(0) + coord_in(1,*) * js(1) + coord_in(2,*) * js(2) z_out = coord_in(0,*) * ks(0) + coord_in(1,*) * ks(1) + coord_in(2,*) * ks(2) endelse return,[x_out, y_out, z_out] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro mkim, coupoff, tmpoff, suz, countf, d2countf, innef, outef, tempf, key_coord = key_coord, lonval = lonval, latval = latval common position, tempsoho, xx_soho, yy_soho, zz_soho, nombre common poterre, xterre_soho, yterre_soho, zterre_soho common rollangle, tfroll, rollval, optroll coupoff = 0. tmpoff = 0. if not keyword_set(key_coord) then begin print, 'Which Coordinate System do you want ?' print, '0 -> SOHO FRAME COORDINATES' print, '1 -> ECLIPTIC COORDINATES' print, '2 -> EQUATORIAL COORDINATES' print, '3 -> INTERSTELLAR WIND FRAME' print, '4 -> GALACTIC COORDINATES' print, '5 -> FIXED SOHO FRAME OF REFERENCE' read, icor endif else icor = fix(key_coord) if (icor lt 0 or icor gt 5) then goto, theend opc = 0 if not keyword_set(lonval) then begin print, ' Longitude: minimum value, maximum value, step size ?' read, lonmin, lonmax, londel endif else begin lonmin = lonval(0) lonmax = lonval(1) londel = lonval(2) endelse if (lonmin lt 0. or lonmax lt 0.) then opc = -1 if not keyword_set(latval) then begin print, ' Latitude: minimum value, maximum value, step size ?' read, latmin, latmax, latdel endif else begin latmin = latval(0) latmax = latval(1) latdel = latval(2) endelse if (londel le 0. or latdel le 0.) then goto, theend nlon = fix((lonmax - lonmin) / londel) if (lonmin + nlon * londel lt lonmax) then nlon = nlon + 1 nlat = fix((latmax - latmin) / latdel) if (latmin + nlat * latdel lt latmax) then nlat = nlat + 1 if (nlon le 0 or nlat le 0) then goto, theend coupoff = fltarr(nlon + 1,nlat,2) coupoff(nlon,0,1) = lonmin coupoff(nlon,1,1) = londel coupoff(nlon,2,1) = latmin coupoff(nlon,3,1) = latdel tmpoff = fltarr(nlon+1,nlat,2) tmpoff(nlon,0,1) = lonmin tmpoff(nlon,1,1) = londel tmpoff(nlon,2,1) = latmin tmpoff(nlon,3,1) = latdel tref = tempf(0) pos_ref = trouve3d(tref,tempsoho,xx_soho,yy_soho,zz_soho) ter_pos = trouve3d(tref,tempsoho,xterre_soho,yterre_soho,zterre_soho) case icor of 0: ter_pos = soho2ecli(pos_ref,ter_pos,-1) 1: ter_pos = ter_pos 2: ter_pos = ecli2equa(ter_pos,1) 3: ter_pos = ecli2wind(ter_pos,1) 4: ter_pos = equa2gala(ecli2equa(ter_pos,1),1) 5: ter_pos = soho2ecli(pos_ref,ter_pos,-1) endcase if optroll eq 1 then begin nrol = n_elements(rollval) if nrol lt 1 then optroll = 0 else begin nrl = where(tfroll le tempf(0),nnl) if nnl gt 0 then nrl = max(nrl) else nrl = 0 endelse endif for ll = 0, n_elements(tempf) - 1 do begin scposit = trouve3d(tempf(ll),tempsoho,xx_soho,yy_soho,zz_soho) out = outef(ll) inn = innef(ll) vout = motoffset(suz,out,inn) out = vout(0) inn = vout(1) tsout = fix(6000 + 40 * out) tsinn = fix(6000 + 40 * inn) if not keeparea(suz,tsinn,tsout) then goto, suivant visee1 = su2soho27(out,inn,suz) anrol = 0. if (icor gt 0 and icor lt 5) then begin if optroll eq 1 then begin suitrol: if nrl lt nrol - 1 then begin if tempf(ll) ge tfroll(nrl + 1) then begin nrl = nrl + 1 goto, suitrol endif endif anrol = rollval(nrl) visee3 = sohoroll2ecli(scposit,visee1,1,anrol) endif else visee3 = soho2ecli(scposit,visee1,1) if icor eq 2 then visee3 = ecli2equa(visee3,1) if icor eq 3 then visee3 = ecli2wind(visee3,1) if icor eq 4 then visee3 = equa2gala(ecli2equa(visee3,1),1) visee1 = visee3 endif if icor eq 5 then begin angll = total(scposit * pos_ref) if angll lt 1. then angll = acos(angll) else angll = 0. vv = visee1 visee1(0,*) = vv(0,*) * cos(angll) - vv(1,*) * sin(angll) visee1(1,*) = vv(0,*) * sin(angll) + vv(1,*) * cos(angll) visee1(2,*) = vv(2,*) endif if opc eq 0 then visee2 = convert3to2deg(visee1) else visee2 = convert3to2deg0(visee1) countjr = 24. * (tempf(ll) - tempf(0)) for rr = 0, 24 do begin vtest = [visee1(0,rr), visee1(1,rr), visee1(2,rr)] vtes1 = ang2vec(vtest,ter_pos) if vtes1 lt 20. then goto, suivantrr jlon = (visee2(0,rr) - lonmin) / londel - 0.5 jlat = (visee2(1,rr) - latmin) / latdel - 0.5 if (jlon ge -0.5 and jlat gt -0.5 and jlon le nlon - 0.5 and jlat lt nlat - 0.5) then begin loni = fix(jlon) dlon = jlon - 1. * loni lons = min([fix(jlon + 1.), nlon - 1]) if jlon lt 0. then begin loni = nlon - 1 lons = 0 dlon = 0. - jlon endif if jlon gt nlon - 1 then begin loni = nlon - 1 lons = 0 dlon = jlon - 1.*loni endif lati = fix(jlat) dlat = jlat - 1. * lati lats = min([fix(jlat + 1.),nlat - 1]) coupoff(loni,lati,0) = coupoff(loni,lati,0) + countf(ll,rr) * (1. - dlon) * (1. - dlat) coupoff(loni,lati,1) = coupoff(loni,lati,1) + d2countf(ll,rr) * (1. - dlon) * (1. - dlat) * (1. - dlon) * (1. - dlat) tmpoff(loni,lati,0) = tmpoff(loni,lati,0) + countjr * (1. - dlon) * (1.-dlat) tmpoff(loni,lati,1) = tmpoff(loni,lati,1) + (1. - dlon) * (1. - dlat) coupoff(loni,lats,0) = coupoff(loni,lats,0) + countf(ll,rr) * (1. - dlon) * dlat coupoff(loni,lats,1) = coupoff(loni,lats,1) + d2countf(ll,rr) * (1. - dlon) * dlat * (1. - dlon) * dlat tmpoff(loni,lats,0) = tmpoff(loni,lats,0) + countjr * (1. - dlon) * dlat tmpoff(loni,lats,1) = tmpoff(loni,lats,1) + (1. - dlon) * dlat coupoff(lons,lati,0) = coupoff(lons,lati,0) + countf(ll,rr) * dlon * (1. - dlat) coupoff(lons,lati,1) = coupoff(lons,lati,1) + d2countf(ll,rr) * dlon * (1. - dlat) * dlon * (1. - dlat) tmpoff(lons,lati,0) = tmpoff(lons,lati,0) + countjr * dlon * (1. - dlat) tmpoff(lons,lati,1) = tmpoff(lons,lati,1) + (1. - dlat) * dlon coupoff(lons,lats,0) = coupoff(lons,lats,0) + countf(ll,rr) * dlon * dlat coupoff(lons,lats,1) = coupoff(lons,lats,1) + d2countf(ll,rr) * dlon * dlat * dlon * dlat tmpoff(lons,lats,0) = tmpoff(lons,lats,0) + countjr * dlon * dlat tmpoff(lons,lats,1) = tmpoff(lons,lats,1) + dlat * dlon endif suivantrr: endfor suivant: endfor for ii = 0, nlat - 1 do begin for jj = 0, nlon - 1 do begin if coupoff(jj,ii,0) gt 0. then begin coupoff(jj,ii,1) = sqrt(coupoff(jj,ii,1)) / coupoff(jj,ii,0) coupoff(jj,ii,0) = coupoff(jj,ii,0) / tmpoff(jj,ii,1) tmpoff(jj,ii,0) = tmpoff(jj,ii,0) / tmpoff(jj,ii,1) endif else begin coupoff(jj,ii,*) = 0. tmpoff(jj,ii,*) = 0. endelse endfor endfor theend: end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function htd2a, imode, su, ht if su eq 1 then p = [1772., 2.36924] else p = [1781., 2.32666] if imode eq 1 then return, 1. * fix(p(0) + p(1) * ht) else return, i = round((ht - p(0)) / p(1) / 10.) * 10 end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro plotimage, coupoff, minicol, maxicol, limites = limites, slon = slon, slat = slat, vco = vco common titredudessin,titreplot info = size(coupoff) if info(0) ne 3 then goto, theend nlon = info(1) - 1 nlat = info(2) if nlat lt 4 then goto, theend lonmin = coupoff(nlon,0,1) londel = coupoff(nlon,1,1) latmin = coupoff(nlon,2,1) latdel = coupoff(nlon,3,1) lonmax = lonmin + nlon * londel latmax = latmin + nlat * latdel if not keyword_set(limites) then begin plonmax = lonmax plonmin = lonmin platmax = latmax platmin = latmin endif else begin plonmax = limites(1) plonmin = limites(0) platmax = limites(3) platmin = limites(2) endelse !p.position = [0.10, 0.05, 0.90, 0.75] if platmin ge 45. then begin soustitre = titreplot + ': North Polar View' map_set, 90., 0,0, /lam, /grid, /noborder, limit = [fix(platmin), -180., 90., 180.] xyouts, 180, 45, '(0,45)', alignment = 0.5 xyouts, 90, 45, ' (90,45)' endif else begin if platmax le -45. then begin soustitre = titreplot + ': South Polar View' map_set, -90., 0,0, /lam, /grid, /noborder, limit = [-90., -180., fix(platmax), 180.] xyouts, 0, -45, '(180,-45)', alignment = 0.5 xyouts, 90, -45, ' (90,-45)' endif else begin if abs(plonmax - plonmin) le 180 then begin soustitre = titreplot map_set, 0.5 * (platmax + platmin), 180. -0.5 * (plonmax + plonmin), 0, /mer, /grid, /noborder, limit = [platmin, 180 - plonmax, platmax, 180 - plonmin] latmoy = 0.5 * (platmax + platmin) xyouts, 180 - plonmax, platmax, string(platmax, format = '(i3)') + ' ', alignment = 1 xyouts, 180 - plonmax, platmin, string(platmin, format = '(i3)') + ' ', alignment = 1 if (platmin lt 0. and platmax gt 0.) then latmoy = 0. xyouts, 180 - plonmax, latmoy, string(plonmax, format = '(i3)') + ' ', alignment = 1 xyouts, 180 - plonmin, latmoy, ' ' + string(plonmin, format = '(i3)') endif else begin soustitre = titreplot map_set, 0., 180. - 0.5 * (plonmax + plonmin), 0, /sin, /grid, /noborder, limit = [platmin, 180 - plonmax, platmax, 180 - plonmin] latmoy = 0.5 * (platmax + platmin) xyouts, 180 - plonmax, platmax, string(platmax, format = '(i3)') + ' ', alignment = 1 xyouts, 180 - plonmax, platmin, string(platmin, format = '(i3)') + ' ', alignment = 1 if (platmin lt 0. and platmax gt 0.) then latmoy = 0. xyouts, 180 - plonmax, latmoy, string(plonmax, format = '(i3)') + ' ', alignment = 1 xyouts, 180 - plonmin, latmoy, ' ' + string(plonmin, format = '(i3)') endelse endelse endelse mincol = float(minicol) maxcol = float(maxicol) colo = (coupoff(*,*,0) - mincol) / (maxcol - mincol) * (!D.table_size - 1) toto = where(colo lt 0., ntoto) if ntoto gt 0 then colo(toto) = 0. toto = where(colo gt (!D.table_size - 1), ntoto) if ntoto gt 0 then colo(toto) = !D.table_size - 1 for il = 0, nlon - 1 do begin for jl = 0, nlat - 1 do begin if coupoff(il,jl,1) gt 0 then begin ipt = 180. - (lonmin + il * londel) igd = 180. - (lonmin + (1 + il) * londel) jpt = latmin + jl * latdel jgd = latmin + (1 + jl) * latdel if ((180 - ipt ge plonmin) and (180 - igd le plonmax)) then if ((jpt ge platmin) and (jgd le platmax)) then polyfill, [ipt, igd, igd, ipt], [jpt, jpt, jgd, jgd], color = colo(il,jl) endif endfor endfor if not keyword_set(vco) then vco=0 if (keyword_set(slon) and keyword_set(slat)) then begin for ll = 0, n_elements(slon) - 1 do begin if ((slon(ll) ge plonmin) and (slon(ll) le plonmax)) then begin if ((slat(ll) ge platmin) and (slat(ll) le platmax)) then begin plots, 180. -slon(ll), slat(ll), psym = 1, symsize = 2, thick = 3, color = (!D.table_size - 1) plots, 180. -slon(ll), slat(ll), psym = 4, symsize = 2, thick = 3, color = 0 endif endif endfor endif if platmin ge 45. then begin map_set, 90, 0, 0, /lam, /grid, /noborder, /noerase, limit = [fix(platmin), -180., 90., 180.] endif else begin if platmax le -45. then begin map_set, -90, 0, 0, /lam, /grid, /noborder, /noerase, limit = [-90., -180., fix(platmax), 180.] endif else begin if abs(plonmax - plonmin) le 180 then begin map_set, 0.5 * (platmax + platmin), 180. - 0.5 * (plonmax + plonmin), 0, /mer, /grid, /noborder, /noerase, limit = [platmin, 180 - plonmax, platmax, 180 - plonmin] endif else begin map_set, 0., 180. - 0.5 * (plonmax + plonmin), 0, /sin, /grid, /noborder, /noerase,limit = [platmin, 180 - plonmax, platmax, 180 - plonmin] endelse endelse endelse !p.position=[0.10,0.88,0.90,0.94] plot, [mincol, maxcol], [0, 1], /nodata, /noerase, xra = [mincol, maxcol], xsty = 1, xticklen = -0.2, xthick = 1, ysty = 4, title = ' Color Scale ', subtitle = soustitre for jj = 0, !D.table_size - 1 do begin xpt = mincol + jj * (maxcol - mincol) / !D.table_size xgd = xpt + (maxcol - mincol) / !D.table_size polyfill, [xpt, xgd, xgd, xpt], [0, 0, 1, 1], color = jj endfor !p.position = [0, 0, 0, 0] theend: end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function analysis, SdState common position, tempsoho, xx_soho, yy_soho, zz_soho, nombre common poterre, xterre_soho, yterre_soho, zterre_soho common rollangle, tfroll, rollval, optroll common vardet, packet, secondes, detsuz, detsu_z, insuz, outsuz, insu_z, $ outsu_z, suzhvps, su_zhvps, suzhcell, su_zhcell, tiz, ti_z, $ outmotz, innmotz, outmot_z, innmot_z, nbre, id, iprim, $ sec_det common header, hearder, ha, hb, hc, hd, he, hf, hg, hh, hi, hj, hk, hl, $ hm, hn, ho, hp, hq, hr, hs, ht common offset, outoffz, innoffz, outoff_z, innoff_z, xout, xinn, nout, ninn common output, tfskz, tfsk_z, cmean, tnord, tsud, flatsuz, flatsu_z, coefns, cal common titredudessin, titreplot status = 0 widget_control, SdState.SdDebug, set_value = '********************************************************************************', /append read_from_fits, SdState.Fitsfile, widget_io = SdState.SdDebug jourdec = (sec_det - 1167609600) / 3600. / 24. nbre = n_elements(secondes) id1 = where(id eq 1 and indgen(nbre) gt 0, nd1) debdec = jourdec(id1(0)) findec = jourdec(max(id1)) widget_control, SdState.SdDebug, set_value = 'Observation period from ' + jourdec2txt(debdec) + ' to '+ jourdec2txt(findec) + '.', /append ;;;;;;;;;;;;;;;;;;;; ; Image definition ; ;;;;;;;;;;;;;;;;;;;; lonmin = 0. lonmax = 360. londel = 1. nlon = 360 latmin = -90. latmax = 90. latdel = 1. nlat = 180 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; SU+Z BaF2 pixel fullsky maps (in ecliptic coordinates) for masking stars (two levels: 10, 12) ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; fsel = file_search(SdState.Savepath,'MasqueBaf2.dat',count = nsel) if nsel eq 0 then begin widget_control, SdState.SdDebug, set_value = 'Cannot find the BaF2 star mask file (MasqueBaf2.dat).', /append widget_control, SdState.SdDebug, set_value = 'Check the IDL save path!', /append goto, theend endif restore, fsel(0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; SoHO orbit and roll attitude ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; fsel = file_search(SdState.Savepath,'orbitsohoidl.dat',count = nsel) if nsel ne 1 then begin widget_control, SdState.SdDebug, set_value = 'Cannot find the SoHO orbit file (orbitsohoidl.dat).', /append widget_control, SdState.SdDebug, set_value = 'Check the IDL save path!', /append goto, theend endif restore, fsel(0) optroll = 0 rollfile = 'rollangle.dat' if (debdec ge 1386. and debdec lt 1418.95) then begin optroll = 1 endif if (debdec ge 1493.76 and debdec lt 1524.97) then begin optroll = 1 endif if debdec ge 3110.54 then begin optroll = 1 endif if optroll eq 1 then begin fsel = file_search(SdState.Savepath,rollfile,count = nsel) if nsel ne 1 then begin widget_control, SdState.SdDebug, set_value = 'Cannot find the SoHO roll file (' + rollfile + ').', /append widget_control, SdState.SdDebug, set_value = 'Check the IDL save path!', /append goto, theend endif widget_control, SdState.SdDebug, set_value = 'SoHO was rolled during this observation.', /append widget_control, SdState.SdDebug, set_value = 'Loading ' + fsel(0) + '...', /append restore, fsel(0) ir = where(tfroll ge debdec and tfroll le findec, nir) if nir ge 1 then begin if ir(0) gt 0 then ir = [ir(0) - 1, ir] if max(ir) lt n_elements(tfroll) - 1 then ir = [ir, max(ir) + 1] rollval = rollval(ir) tfroll = tfroll(ir) endif else begin b1 = where(tfroll le debdec, nb1) b2 = where(tfroll gt findec, nb2) if nb1 gt 0 then begin if nb2 gt 0 then ir = [max(b1), min(b2)] else ir = [max(b1)] rollval = rollval(ir) tfroll = tfroll(ir) endif endelse endif ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Load the corrections on motors positions ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; fsel = file_search(SdState.Savepath,'offset+Zidl.dat',count = nsel) if nsel ne 1 then begin widget_control, SdState.SdDebug, set_value = 'Cannot find the offset+Zidl.dat file.', /append widget_control, SdState.SdDebug, set_value = 'Check the IDL save path!', /append goto, theend endif restore, fsel(0) fsel = file_search(SdState.Savepath,'offset-Zidl.dat',count = nsel) if nsel ne 1 then begin widget_control, SdState.SdDebug, set_value = 'Cannot find the offset-Zidl.dat file.', /append widget_control, SdState.SdDebug, set_value = 'Check the IDL save path!', /append goto, theend endif restore, fsel(0) ;;;;;;;;;;;;;;;;;;;; ; Graphical output ; ;;;;;;;;;;;;;;;;;;;; widget_control, SdState.SdMotor1, get_value = wdmc1 widget_control, SdState.SdMotor2, get_value = wdmc2 widget_control, SdState.SdFlatfield1, get_value = wdflfi1 widget_control, SdState.SdFlatfield2, get_value = wdflfi2 widget_control, SdState.SdNsc, get_value = wdnsc widget_control, SdState.SdImage, get_value = wdimg ;;;;;;;;;;;;;;;;;;;;;; ; SU+Z data analysis ; ;;;;;;;;;;;;;;;;;;;;;; widget_control, SdState.SdDebug, set_value = '********', /append widget_control, SdState.SdDebug, set_value = 'SU+Z', /append if n_elements(outsuz) gt 1 then begin jourdec1 = 0 cleanid, 1, id, jourdec, jourdec1, outsuz, outmotz, insuz, innmotz, detsuz, suzhcell, suzhvps, tiz, widget_io = SdState.SdDebug endif if n_elements(outsuz) gt 1 then begin widget_control, SdState.SdDebug, set_value = 'Found ' + string(n_elements(outsuz), format = '(i5)') + ' packets.', /append wset, wdmc1 !p.multi = [0, 1, 2] avefixmot, 1, 'Outer', outsuz, outmotz, /plot avefixmot, 1, 'Inner', insuz, innmotz, /plot tfskz = jourdec1(0) ; Correct the noisy pixel 16: noisepix = 39.2181 - 0.00321215 * tfskz detsuz(*,16) = detsuz(*,16) - noisepix widget_control, SdState.SdDebug, set_value = 'Noise pixel 16: ' + string(noisepix, format = '(f6.2)') + ' cps' , /append suzhcell = suzhcell AND '38'X if (n_elements(where(suzhvps gt 2000. and suzhcell ge 8)) gt 1) then begin widget_control, SdState.SdDebug, set_value = 'SU+Z Hcell active.', /append goto, theend endif ; Check the pixel 26 for particle event: widget_control, SdState.SdDebug, set_value = 'Pixel 26, median counting rate: ' + string(median(detsuz(*,26)), format = '(f6.2)') + ' cps' , /append pk = where (abs(outmotz - 6000) le 4000 and abs(innmotz - 6000) le 4000 and suzhvps gt 1900. and suzhcell eq 0,npk) if npk gt 0 then begin countf = fltarr(npk,27) d2countf = fltarr(npk,27) for k = 0, npk - 1 do begin for ii = 0, 26 do begin countf(k,ii) = float(detsuz(pk(k),ii)) d2countf(k,ii) = countf(k,ii) / tiz(pk(k)) endfor endfor innef = (insuz(pk) - 6000.) / 40. outef = (outsuz(pk) - 6000.) / 40. tempf = jourdec1(pk) wset, wdflfi1 flatsuz = flficomp(1,countf,d2countf,innef,outef,tempf,masque10) widget_control, SdState.SdDebug, set_value = 'Flatfield corrections:', /append widget_control, SdState.SdDebug, set_value = 'pixel ratio sigma chi2 / n.d.f.', /append for k = 0, 24 do if k ne 12 then $ widget_control, SdState.SdDebug, set_value = $ ' ' + string(k, format = '(i02)') + ' ' + string(flatsuz(k,0), format = '(f4.2)') + ' ' + string(flatsuz(k,1), format = '(f4.2)') + ' ' + string(flatsuz(k,2), format = '(f4.2)'), /append for k = 0, npk - 1 do begin for ii = 0, 24 do begin countf(k,ii) = float(detsuz(pk(k),ii)) * flatsuz(ii,0) d2countf(k,ii) = countf(k,ii)^2 * ((flatsuz(ii,1) / flatsuz(ii,0))^2 + (1. / (float(detsuz(pk(k),ii)) * tiz(pk(k))))) endfor endfor mkim, cnord, tnord, 1, countf, d2countf, innef, outef, tempf, key_coord = 1, lonval = [0., 360., 1.], latval = [-90., 90., 1.] endif endif else widget_control, SdState.SdDebug, set_value = 'No data for SU+Z', /append ;;;;;;;;;;;;;;;;;;;;;; ; SU-Z data analysis ; ;;;;;;;;;;;;;;;;;;;;;; widget_control, SdState.SdDebug, set_value = '********', /append widget_control, SdState.SdDebug, set_value = 'SU-Z', /append if n_elements(outsu_z) gt 1 then begin jourdec2 = 0 cleanid, -1, id, jourdec, jourdec2, outsu_z, outmot_z, insu_z, innmot_z, detsu_z, su_zhcell, su_zhvps, ti_z, widget_io = SdState.SdDebug endif if n_elements(outsu_z) gt 1 then begin widget_control, SdState.SdDebug, set_value = 'Found ' + string(n_elements(outsu_z), format = '(i5)') + ' packets.', /append checkoutinn, outsu_z, outmot_z, insu_z, innmot_z, widget_io = SdState.SdDebug wset, wdmc2 !p.multi = [0, 1, 2] avefixmot, -1, 'Outer', outsu_z, outmot_z, /plot avefixmot, -1, 'Inner', insu_z, innmot_z, /plot tfsk_z = jourdec2(0) su_zhcell = su_zhcell AND '38'X if (n_elements(where(su_zhvps gt 2000. and su_zhcell ge 8)) gt 1) then begin widget_control, SdState.SdDebug, set_value = 'SU-Z Hcell active.', /append goto, theend endif pk = where (abs(outmot_z - 6000) le 4000 and abs(innmot_z - 6000) le 4000 and su_zhvps gt 1900. and su_zhcell eq 0,npk) if npk gt 0 then begin countf = fltarr(npk,27) d2countf = fltarr(npk,27) for k = 0, npk - 1 do begin for ii = 0, 26 do begin countf(k,ii) = float(detsu_z(pk(k),ii)) d2countf(k,ii) = countf(k,ii) / ti_z(pk(k)) endfor endfor innef = (insu_z(pk) - 6000.) / 40. outef = (outsu_z(pk) - 6000.) / 40. innec = innef * 0. outec = outef * 0. b1 = where(outef lt 0.,nb1) if nb1 gt 0 then begin innec(b1) = -4.11 - 0.0497 * outef(b1) - 0.000321 * outef(b1)^2 outec(b1) = -0.181 - 0.0137 * outef(b1) innef(b1) = innef(b1) + innec(b1) outef(b1) = outef(b1) + outec(b1) endif b2 = where(outef gt 0.,nb2) if nb2 gt 0 then begin innec(b2) = 0.17 - 0.0484 * outef(b2) + 0.000312 * outef(b2)^2 outec(b2) = -0.125-0.0103 * outef(b2) innef(b2) = innef(b2) + innec(b2) outef(b2) = outef(b2) + outec(b2) endif tempf = jourdec2(pk) wset, wdflfi2 flatsu_z = flficomp(-1,countf,d2countf,innef,outef,tempf,masque10) widget_control, SdState.SdDebug, set_value = 'Flatfield corrections:', /append widget_control, SdState.SdDebug, set_value = 'pixel ratio sigma chi2 / n.d.f.', /append for k = 0, 24 do if k ne 12 then $ widget_control, SdState.SdDebug, set_value = $ ' ' + string(k, format = '(i02)') + ' ' + string(flatsu_z(k,0), format = '(f4.2)') + ' ' + string(flatsu_z(k,1), format = '(f4.2)') + ' ' + string(flatsu_z(k,2), format = '(f4.2)'), /append for k = 0, npk - 1 do begin for ii = 0, 24 do begin countf(k,ii) = float(detsu_z(pk(k),ii)) * flatsu_z(ii,0) d2countf(k,ii) = countf(k,ii)^2 * ((flatsu_z(ii,1) / flatsu_z(ii,0))^2 + (1. / (float(detsu_z(pk(k),ii)) * ti_z(pk(k))))^2) endfor endfor mkim, csud, tsud, -1, countf, d2countf, innef, outef, tempf, key_coord = 1, lonval = [0., 360., 1.], latval = [-90., 90., 1.] endif endif else widget_control, SdState.SdDebug, set_value = 'No data for SU-Z', /append ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; SU+Z absolute calibration ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; widget_control, SdState.SdDebug, set_value = '********', /append widget_control, SdState.SdDebug, set_value = 'Setting absolute calibration for SU+Z:', /append cal = [-1., 0.] if debdec lt 2306. then dhvr1 = 110 else dhvr1 = 120 if htd2a(0,1,mean(suzhvps)) eq dhvr1 then begin if debdec le 1269 then cal = [1.190, 0.120] if debdec ge 1388 and debdec le 2147 then cal = [3.037, 0.374] if debdec ge 2160 and debdec le 2306 then cal = [3.924, 0.551] if debdec ge 2307 then cal = [2.512, 0.403] endif else widget_control, SdState.SdDebug, set_value = 'SU+Z high voltage not at its nominal value.', /append if cal(0) lt 0 then begin widget_control, SdState.SdDebug, set_value = 'Data are in count by second, no calibration applied !', /append cal(0) = 1. endif else widget_control, SdState.SdDebug, set_value = 'Overall calibration factor: ' + string(cal(0), format = '(f4.2)') + ' +/- ' + string(cal(1), format = '(f4.2)') + ' R/cps.', /append ;;;;;;;;;;;;;;;;;;;;;;;;;; ; North/South coeficient ; ;;;;;;;;;;;;;;;;;;;;;;;;;; cn = fltarr(nlon,nlat) en = fltarr(nlon,nlat) cs = fltarr(nlon,nlat) es = fltarr(nlon,nlat) for ilon = 0, nlon - 1 do begin for ilat = 0, nlat - 1 do begin if cnord(ilon,ilat,1) gt 0. then begin cn(ilon,ilat) = cnord(ilon,ilat,0) * cal(0) * masque10(ilon,ilat) en(ilon,ilat) = sqrt((cal(1) / cal(0)) ^ 2 + cnord(ilon,ilat,1) ^ 2) * masque10(ilon,ilat) endif if csud(ilon,ilat,1) gt 0. then begin cs(ilon,ilat) = csud(ilon,ilat,0) * cal(0) * masque10(ilon,ilat) es(ilon,ilat) = sqrt((cal(1) / cal(0)) ^ 2 + csud(ilon,ilat,1) ^ 2) * masque10(ilon,ilat) endif endfor endfor com = where(en gt 0. and es gt 0. and cn lt 1000., ncom) coefns = (ncom gt 0 ? fpr(cs(com),es(com),cn(com),en(com),1.) : -1) wset, wdnsc !p.multi = 0 plot, cn(com), cs(com), psy = 3, xtitle = 'North [cps]', ytitle = 'South [cps]', xrange = [min(cn(com)), max(cn(com))], yrange = [min(cs(com)), max(cs(com))] pts = [min(cn(com)), max(cn(com))] oplot, pts, pts / coefns(0) if n_elements(coefns) ne 3 then begin widget_control, SdState.SdDebug, set_value = 'Cannot calculate the N/S coefficient!', /append coefns = [1., 0] endif else widget_control, SdState.SdDebug, set_value = 'N/S coefficient: ' + string(coefns(0), format = '(f3.1)') + ' +/- ' + string(coefns(1), format = '(f3.1)') + '.', /append ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Gather both half maps together ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; cmean = cnord for ilon = 0, nlon - 1 do begin for ilat = 0, nlat - 1 do begin if cnord(ilon,ilat,1) gt 0. then begin cnord(ilon,ilat,0) = cnord(ilon,ilat,0) * cal(0) cnord(ilon,ilat,1) = sqrt((cal(1) / cal(0)) ^ 2 + cnord(ilon,ilat,1) ^ 2) endif if csud(ilon,ilat,1) gt 0. then begin csud(ilon,ilat,0) = csud(ilon,ilat,0) * cal(0) * coefns(0) csud(ilon,ilat,1) = sqrt((cal(1) / cal(0)) ^ 2 + csud(ilon,ilat,1) ^ 2 + (coefns(1) / coefns(0)) ^ 2) endif if (cnord(ilon,ilat,1) gt 0. and csud(ilon,ilat,1) gt 0.) then begin cmean(ilon,ilat,0) = 0.5 * (cnord(ilon,ilat,0) + csud(ilon,ilat,0)) cmean(ilon,ilat,1) = sqrt((cnord(ilon,ilat,1) * cnord(ilon,ilat,0)) ^2 + (csud(ilon,ilat,1) * csud(ilon,ilat,0)) ^ 2) / (cnord(ilon,ilat,0) + csud(ilon,ilat,0)) endif if (cnord(ilon,ilat,1) gt 0. and csud(ilon,ilat,1) le 0.) then cmean(ilon,ilat,0:1) = cnord(ilon,ilat,0:1) if (cnord(ilon,ilat,1) le 0. and csud(ilon,ilat,1) gt 0.) then cmean(ilon,ilat,0:1) = csud(ilon,ilat,0:1) endfor endfor titreplot = jourdec2txt(debdec) wset, wdimg plotimage, cmean, 0., 1000. widget_control, SdState.SdDebug, set_value = 'Done.', /append status = 1 theend: return, status end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro SdEvent, Event common output, tfskz, tfsk_z, cmean, tnord, tsud, flatsuz, flatsu_z, coefns, cal common titredudessin, titreplot widget_control, get_uvalue = SdState, Event.Top, /no_copy widget_control, SdState.SdFitspath, get_value = dummy SdState.FitsPath = dummy widget_control, SdState.SdIdlsavepath, get_value = dummy SdState.Savepath = dummy case Event.Id of SdState.SdOpen: begin SdState.Fitsfile = Dialog_Pickfile(/read, path = SdState.FitsPath, filter = 'fsky??????.fits', dialog_parent = Event.Top, title = 'Please select a fits data file') widget_control, SdState.SdBase, tlb_set_title = 'SoHO/SWAN fullsky display - ' + file_basename(SdState.Fitsfile) widget_control, SdState.SdBaseMotor1, tlb_set_title = 'SU+Z motors correction - ' + file_basename(SdState.Fitsfile) widget_control, SdState.SdBaseMotor2, tlb_set_title = 'SU-Z motors correction - ' + file_basename(SdState.Fitsfile) widget_control, SdState.SdBaseFlatfield1, tlb_set_title = 'SU+Z flatfield corrections - ' + file_basename(SdState.Fitsfile) widget_control, SdState.SdBaseFlatfield2, tlb_set_title = 'SU-Z flatfield corrections - ' + file_basename(SdState.Fitsfile) widget_control, SdState.SdBaseNsc, tlb_set_title = 'N/S coef. - ' + file_basename(SdState.Fitsfile) widget_control, SdState.SdBaseImage, tlb_set_title = 'Image display - ' + file_basename(SdState.Fitsfile) endcase SdState.SdSavebin: begin fileout = dialog_pickfile(file = file_basename(SdState.Fitsfile, '.fits'), /write, default_extension = 'bin') save, tfskz, tfsk_z, cmean, tnord, tsud, titreplot, flatsuz, flatsu_z, coefns, cal, filename = fileout widget_control, SdState.SdDebug, set_value = '********', /append widget_control, SdState.SdDebug, set_value = 'Output saved in ' + fileout + '.', /append endcase SdState.SdSavegif: begin fileout = dialog_pickfile(file = file_basename(SdState.Fitsfile, '.fits'), /write, default_extension = 'gif') widget_control, SdState.SdImage, get_value = wdimg dname = !d.name set_plot, 'z' plotimage, cmean, 0., 1000. tvlct, r, g, b, /get write_gif, fileout, tvrd(), r, g, b set_plot, dname widget_control, SdState.SdDebug, set_value = '********', /append widget_control, SdState.SdDebug, set_value = 'Image saved in ' + fileout + '.', /append endcase SdState.SdVersion: begin widget_control, SdState.SdDebug, set_value = '********', /append widget_control, SdState.SdDebug, set_value = 'This is swand $Revision: 1.9 $', /append endcase SdState.SdRun: SdState.Done = analysis(SdState) SdState.SdMenuMotor1: widget_control, SdState.SdBaseMotor1, map = (Widget_Info(SdState.SdBaseMotor1,/map) ? 0 : 1) SdState.SdMenuMotor2: widget_control, SdState.SdBaseMotor2, map = (Widget_Info(SdState.SdBaseMotor2,/map) ? 0 : 1) SdState.SdMenuFlatfield1: widget_control, SdState.SdBaseFlatfield1, map = (Widget_Info(SdState.SdBaseFlatfield1,/map) ? 0 : 1) SdState.SdMenuFlatfield2: widget_control, SdState.SdBaseFlatfield2, map = (Widget_Info(SdState.SdBaseFlatfield2,/map) ? 0 : 1) SdState.SdMenuNsc: widget_control, SdState.SdBaseNsc, map = (Widget_Info(SdState.SdBaseNsc,/map) ? 0 : 1) SdState.SdMenuImage: widget_control, SdState.SdBaseImage, map = (Widget_Info(SdState.SdBaseImage,/map) ? 0 : 1) SdState.SdSetpaths: widget_control, SdState.SdBaseSetpaths, map = (Widget_Info(SdState.SdBaseSetpaths,/map) ? 0 : 1) else: endcase widget_control, SdState.SdRun, sensitive = (SdState.Fitsfile eq '' ? 0 : 1) widget_control, SdState.SdSavebin, sensitive = (SdState.Done eq '' ? 0 : 1) widget_control, SdState.SdSavegif, sensitive = (SdState.Done eq '' ? 0 : 1) if Event.Id eq SdState.SdExit then begin widget_control, SdState.SdBase, /destroy widget_control, SdState.SdBaseMotor1, /destroy widget_control, SdState.SdBaseMotor2, /destroy widget_control, SdState.SdBaseFlatfield1, /destroy widget_control, SdState.SdBaseFlatfield2, /destroy widget_control, SdState.SdBaseNsc, /destroy widget_control, SdState.SdBaseImage, /destroy widget_control, SdState.SdBaseSetpaths, /destroy endif else widget_control, set_uvalue = SdState, Event.Top, /no_copy end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro swand compile_opt hidden on_error, 2 if xregistered("swand") ne 0 then return case 1 of !d.name eq 'WIN' : device, decomposed = 0, retain = 2 !d.name eq 'X' or !d.name eq 'MAC' : device, true_color = 24, decomposed = 0, retain = 2 endcase loadct, 38 defaultfitspath = '' defaultSavepath = '' SdBase = widget_base(group_leader = wGroup, uname = 'SdBase', title = 'SoHO/SWAN fullsky display', xoffset = 10, yoffset = 10, space = 3, xpad = 3, ypad = 3, tlb_frame_attr = 1, mbar = SdMenu) SdBaseMotor1 = widget_base(group_leader = wGroup, uname = 'SdBaseMotor1', title = 'SU+Z motors correction', map = 0, tlb_frame_attr = 8) SdBaseMotor2 = widget_base(group_leader = wGroup, uname = 'SdBaseMotor2', title = 'SU-Z motors correction', map = 0, tlb_frame_attr = 8) SdBaseFlatfield1 = widget_base(group_leader = wGroup, uname = 'SdBaseFlatfield1', title = 'SU+Z flatfield corrections', map = 0, tlb_frame_attr = 8) SdBaseFlatfield2 = widget_base(group_leader = wGroup, uname = 'SdBaseFlatfield2', title = 'SU-Z flatfield corrections', map = 0, tlb_frame_attr = 8) SdBaseNsc = widget_base(group_leader = wGroup, uname = 'SdBaseNsc', title = 'N/S coef.', map = 0, tlb_frame_attr = 8) SdBaseImage = widget_base(group_leader = wGroup, uname = 'SdBaseImage', title = 'Image display', map = 1, tlb_frame_attr = 8, xoffset = 400, yoffset = 200) SdBaseSetpaths = widget_base(group_leader = wGroup, uname = 'SdBaseSetpaths', xoffset = 80, yoffset = 440, title = 'Search paths', map = 0, tlb_frame_attr = 8) SdMenuFile = widget_button(SdMenu, uname = 'SdMenuFile', /menu, value = 'File') SdOpen = widget_button(SdMenuFile, uname = 'SdOpen' , value = 'Open', accelerator = "Ctrl+O") SdSetpaths = widget_button(SdMenuFile, uname = 'SdSetpaths', value = 'Set paths...', /separator, accelerator = "Ctrl+P") SdRun = widget_button(SdMenuFile, uname = 'SdRun', value = 'Proceed', sensitive = 0, /separator, accelerator = "Ctrl+R") SdSavebin = widget_button(SdMenuFile, uname = 'SdSavebin' , value = 'Save bin', sensitive = 0, /separator, accelerator = "Ctrl+S") SdSavegif = widget_button(SdMenuFile, uname = 'SdSavegif' , value = 'Save gif', sensitive = 0, accelerator = "Ctrl+G") SdExit = widget_button(SdMenuFile, uname = 'SdExit' , value = 'Exit', /separator, accelerator = "Ctrl+Q") SdMenuWindow = widget_button(SdMenu, uname = 'SdMenuWindow', /menu, value = 'Window') SdMenuSU1 = widget_button(SdMenuWindow, uname = 'SdMenuSU1', value = 'SU+Z', /menu) SdMenuSU2 = widget_button(SdMenuWindow, uname = 'SdMenuSU2', value = 'SU-Z', /menu) SdMenuMotor1 = widget_button(SdMenuSU1, uname = 'SdMenuMotor1', value = 'Motors correction') SdMenuMotor2 = widget_button(SdMenuSU2, uname = 'SdMenuMotor2', value = 'Motors correction') SdMenuFlatfield1 = widget_button(SdMenuSU1, uname = 'SdMenuFlatfield1', value = 'Flatfield corrections') SdMenuFlatfield2 = widget_button(SdMenuSU2, uname = 'SdMenuFlatfield2', value = 'Flatfield corrections') SdMenuNsc = widget_button(SdMenuWindow, uname = 'SdMenuSnc', value = 'N/S coef.') SdMenuImage = widget_button(SdMenuWindow, uname = 'SdMenuImage', value = 'Image display') SdMenuHelp = widget_button(SdMenu, uname = 'SdMenuHelp', /menu, /help, value ='Help') SdHowTo = widget_button(SdMenuHelp, uname = 'SdHowTo', sensitive = 0, value = 'HowTo...') SdVersion = widget_button(SdMenuHelp, uname = 'SdVersion', /separator, value = 'Version') dummy = widget_label(SdBaseSetpaths, value='fits data:') SdFitspath = widget_text(SdBaseSetpaths, uname = 'SdFitspath', xoffset = 100, /editable, xsize = 40, value = defaultfitspath) dummy = widget_label(SdBaseSetpaths, yoffset = 40, value = 'IDL save files:') SdIdlsavepath = widget_text(SdBaseSetpaths, uname = 'SdIdlsave', xoffset = 100, yoffset = 40, /editable, xsize = 40, value = defaultSavepath) SdDebug = widget_text(SdBase, uname = 'SdDebug', xsize = 80, ysize = 25, /scroll) SdMotor1 = Widget_Draw(SdBaseMotor1, uname = 'SdMotor1', scr_xsize = 500, scr_ysize = 500) SdMotor2 = Widget_Draw(SdBaseMotor2, uname = 'SdMotor2', scr_xsize = 500, scr_ysize = 500) SdFlatfield1 = Widget_Draw(SdBaseFlatfield1, uname = 'SdFlatfield1', scr_xsize = 1000, scr_ysize = 800) SdFlatfield2 = Widget_Draw(SdBaseFlatfield2, uname = 'SdFlatfield2', scr_xsize = 1000, scr_ysize = 800) SdNsc = Widget_Draw(SdBaseNsc, uname = 'SdNsc', scr_xsize = 300, scr_ysize = 300) SdImage = Widget_Draw(SdBaseImage, uname = 'SdImage', scr_xsize = 800, scr_ysize = 600) widget_control, SdBase, /realize widget_control, SdBaseMotor1, /realize widget_control, SdBaseMotor2, /realize widget_control, SdBaseFlatfield1, /realize widget_control, SdBaseFlatfield2, /realize widget_control, SdBaseNsc, /realize widget_control, SdBaseImage, /realize widget_control, SdBaseSetpaths, /realize SdState = { $ ; Default value: Fitsfile : '', $ FitsPath : defaultfitspath, $ Savepath : defaultSavepath, $ Done : 0, $ ; Main window: SdBase : SdBase, $ SdDebug : SdDebug, $ ; SU+Z motor corrections window: SdMenuMotor1 : SdMenuMotor1, $ SdBaseMotor1 : SdBaseMotor1, $ SdMotor1 : SdMotor1, $ ; SU+Z flatfield corrections fit: SdMenuFlatfield1 : SdMenuFlatfield1, $ SdBaseFlatfield1 : SdBaseFlatfield1, $ SdFlatfield1 : SdFlatfield1, $ ; SU-Z motor corrections window: SdMenuMotor2 : SdMenuMotor2, $ SdBaseMotor2 : SdBaseMotor2, $ SdMotor2 : SdMotor2, $ ; SU-Z flatfield corrections fit: SdMenuFlatfield2 : SdMenuFlatfield2, $ SdBaseFlatfield2 : SdBaseFlatfield2, $ SdFlatfield2 : SdFlatfield2, $ ; N/S coeficient fit: SdMenuNsc : SdMenuNsc, $ SdBaseNsc : SdBaseNsc, $ SdNsc : SdNsc, $ ; Image window: SdMenuImage : SdMenuImage, $ SdBaseImage : SdBaseImage, $ SdImage : SdImage, $ ; Set paths window: SdBaseSetpaths : SdBaseSetpaths, $ SdFitspath : SdFitspath, $ SdIdlsavepath : SdIdlsavepath, $ ; File Menu: SdOpen : SdOpen, $ SdRun : SdRun, $ SdSetpaths : SdSetpaths, $ SdSavebin : SdSavebin, $ SdSavegif : SdSavegif, $ SdExit : SdExit, $ ; Help Menu: SdHowTo : SdHowTo, $ SdVersion : SdVersion $ } widget_control, SdBase, set_uvalue = SdState xmanager, "swand", SdBase, event_handler = "SdEvent", group_leader = wGroup, /no_block end