-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathscc_measure3.pro
1474 lines (1451 loc) · 51.5 KB
/
scc_measure3.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
;+
; $Id: scc_measure3.pro, 2024/04/24 sahade Exp $
; DOI 10.5281/zenodo.13951841
; Project : STEREO - SECCHI
;
; Name : SCC_MEASURE3
;
; Purpose : Measure 3D coordinates from three STEREO images
;
; Category : STEREO, SECCHI, Coordinates
;
; Explanation : Widget application to allow a user to select with the cursor a
; feature appearing in three images. The selected coordinates are
; used to locate the point in three-dimensional space.
;
; The user first selects a point on one image. The program then
; displays a line on the other two images representing the
; line-of-sight from the first image. The user then selects the
; point along this line with the same feature, in any of the other images.
; The 3D coordinates are then calculated and displayed within the widget
; as (Earth-based) Stonyhurst heliographic longitude and
; latitude, along with radial distance in solar radii.
;
; In essence, this is a demonstration procedure showing how to
; derive 3D information from multi-viewpoint data.
; Tested with: EUVI A & B, EIT, AIA and EUI
;
; Syntax : SCC_MEASURE3, IMAGE_AHEAD, IMAGE_BEHIND, IMAGE_OTHER, INDEX_AHEAD, INDEX_BEHIND, INDEX_OTHER
;
;
; or
;
; SCC_MEASURE, FILE_AHEAD, FILE_BEHIND, FILE_OTHER
; Only works if Other can be read by sccreadfits
;
; Inputs : IMAGE_AHEAD = Image for the Ahead observatory
; IMAGE_BEHIND = Image for the Behind observatory
; IMAGE_OTHER = Image for the Other observatory
; INDEX_AHEAD = FITS index structure for Ahead
; INDEX_BEHIND = FITS index structure for Behind
; INDEX_OTHER = FITS index structure for Other
;
; Prefered method
;
; Not tested:
; Opt. Inputs : FILE_AHEAD = Filename for the Ahead observatory.
; FILE_BEHIND = Filename for the Behind observatory
; FILE_OTHER = Filename for other observatory
;
; ;
;
; Outputs : None.
;
; Opt. Outputs: None.
;
; Keywords : WSIZE = Window size to use for displaying the images. The
; default is 512.
;
; OUTFILE = Name of file to store results in. Use the "Store"
; button to write each selected measurement into this
; file. The stored points will also be displayed on
; the images. Use the "Clear stored" button to clear
; these points out from memory. At the same time, a
; blank line will be written to the output file to mark
; the separation between collections of points.
;
; APPEND = Append to existing output file.
;
; FORCESAVE = If set, saves last results in given output
; file upon exit, even if manual 'Save' has
; or has not been used. Has no effect if no
; output file was defined.
;
; CROP = If set, restricts chosing of points to the
; inner region so that a crop to half size
; can later be done (outside of scc_measure)
;
; NO_BLOCK = Call XMANAGER with /NO_BLOCK
;
; /DEBUG print diagnostic information
;
;
; Calls : SCCREADFITS, FITSHEAD2WCS, EXPTV, SCC_MEASURE3_EVENT,
; SCC_MEASURE3_SELECT_1ST, SCC_MEASURE3_SELECT_2ND, BOOST_ARRAY,
; SCC_MEASURE3_REDISPLAY, SSC_MEASURE3_REPLOT, SSC_MEASURE3_CLEANUP,
; SETWINDOW, GET_TV_SCALE, WCS_GET_COORD, CONVERT_SUNSPICE_COORD
;
; Common : SCC_MEASURE3 = Internal common block used by the widget program.
;
; Restrictions: None.
;
; Side effects: None.
;
;
; Adapted from SCC_MEASURE (William Thompson)
;
; Written : Version 1, 26-Jan-2023, Abril Sahade
;-
;
;==============================================================================
;
pro scc_measure3_redisplay, window=k_window
;
; Procure to redisplay the windows after various operations such as changing
; the color table.
;
common scc_measure3, ahead, h_ahead, behind, h_behind, eother, h_eother, $
wcs_ahead, wcs_behind, wcs_eother, win_left, win_right, win_other, $
image_left, image_right, image_other, lon_wid, lat_wid, rsun_wid, $
zoom_wid, unzoom_wid, maxz, pix_paraml, pix_paramr, pix_paramo, $
pix_left, pix_right, pix_other, heeq_left, heeq_right, heeq_other, $
in_progress, win_last, lzoom, rzoom, ozoom, lfrz_wid, rfrz_wid, $
ofrz_wid, color, color1, color2, lmin, lmax, rmin, rmax, omin, omax, $
lmin_wid, lmax_wid, rmin_wid, rmax_wid, omin_wid, omax_wid, $
subimage_left, subimage_right, subimage_other, $
origin_left, origin_right, origin_other, outlun, out_wid, $
clear_wid, oplot_wid, n_stored, store_left, store_right, store_other, $
pixpair, pixforce, roi, debugon
;
if n_elements(k_window) eq 1 then window = k_window else window = -1
;
; If the left (or no) window was selected, then redisplay the left window.
;
if (window eq win_left) or (window eq -1) then begin
setwindow, win_left
get_tv_scale, sx, sy, mx, my, jx, jy
tv, image_left, jx, jy
if n_stored gt 0 then $
plots, store_left[0,*], store_left[1,*], psym=1, color=color
endif
;
; If the right (or no) window was selected, then redisplay the right window.
;
if (window eq win_right) or (window eq -1) then begin
setwindow, win_right
get_tv_scale, sx, sy, mx, my, jx, jy
tv, image_right, jx, jy
if n_stored gt 0 then $
plots, store_right[0,*], store_right[1,*], psym=1, color=color
endif
;
;; If the other (or no) window was selected, then redisplay the right window.
;
if (window eq win_other) or (window eq -1) then begin
setwindow, win_other
get_tv_scale, sx, sy, mx, my, jx, jy
tv, image_other, jx, jy
if n_stored gt 0 then $
plots, store_other[0,*], store_other[1,*], psym=1, color=color
endif
;
; If no particular window was selected, and some points have been selected,
; then recreate the plots on top of the image.
;
if (window eq -1) and (win_last ne '') then begin
;
; If the point selection is in progress, then plot the point on the 1st
; window, and the line on the 2nd and other window.
;
if in_progress then begin
if win_last eq 'LEFT' then begin
setwindow, win_left
plots, pix_left[0], pix_left[1], psym=1, symsize=3, color=color
setwindow, win_right
plots, pix_paramr[0,*], pix_paramr[1,*], color=color
x = [0,wcs_behind.naxis[1]-1]
paramr = poly_fit(pix_paramr[0,*], pix_paramr[1,*], 1)
plots, x, poly(x, paramr), line=1, color=color
setwindow, win_other
plots, pix_paramo[0,*], pix_paramo[1,*], color=color
x2 = [0,wcs_eother.naxis[1]-1]
paramo = poly_fit(pix_paramo[0,*], pix_paramo[1,*], 1)
plots, x2, poly(x2, paramo), line=1, color=color
end else if win_last eq 'RIGHT' then begin
setwindow, win_right
plots, pix_right[0], pix_right[1], psym=1, symsize=3, color=color
setwindow, win_left
plots, pix_paraml[0,*], pix_paraml[1,*], color=color
x = [0,wcs_ahead.naxis[1]-1]
paraml = poly_fit(pix_paraml[0,*], pix_paraml[1,*], 1)
plots, x, poly(x, paraml), line=1, color=color
setwindow, win_other
plots, pix_paramo[0,*], pix_paramo[1,*], color=color
x2 = [0,wcs_eother.naxis[1]-1]
paramo = poly_fit(pix_paramo[0,*], pix_paramo[1,*], 1)
plots, x2, poly(x2, paramo), line=1, color=color
end else begin
setwindow, win_other
plots, pix_other[0], pix_other[1], psym=1, symsize=3, color=color
setwindow, win_left
plots, pix_paraml[0,*], pix_paraml[1,*], color=color
x = [0,wcs_ahead.naxis[1]-1]
paraml = poly_fit(pix_paraml[0,*], pix_paraml[1,*], 1)
plots, x, poly(x, paraml), line=1, color=color
setwindow, win_right
plots, pix_paramr[0,*], pix_paramr[1,*], color=color
x2 = [0,wcs_behind.naxis[1]-1]
paramr = poly_fit(pix_paramr[0,*], pix_paramr[1,*], 1)
plots, x2, poly(x2, paramr), line=1, color=color
endelse
;
; Otherwise, plot both points.
;
end else begin
setwindow, win_left
plots, pix_left[0], pix_left[1], psym=1, symsize=3, color=color
setwindow, win_right
plots, pix_right[0], pix_right[1], psym=1, symsize=3, color=color
setwindow, win_other
;x = [0,wcs_eother.naxis[1]-1]
;plots, pix_param2[0,*], pix_param2[1,*], color=color
;plots, x, poly(x, param2), line=1, color=color
;plots, pix_param3[0,*], pix_param3[1,*], color=color1
;plots, x, poly(x, param3), line=1, color=color1
plots, pix_other[0], pix_other[1], psym=1, symsize=3, color=color ;
endelse
endif
;
end
;
;==============================================================================
;
pro scc_measure3_oplot_from_file ;without modifications from 2img version
;
; Procedure to overplot traces already saved in the output file.
;
common scc_measure3
;
; Read in the file.
;
infile = (fstat(outlun)).name
openr, inlun, infile, /get_lun
delvarx, table, iseg
line = 'string'
noffset0 = 0
values = dblarr(7)
while not eof(inlun) do begin
readf, inlun, line
;
; Only use lines containing 7 numbers.
;
on_ioerror, finish_segment
reads, line, values
boost_array, table, values
goto, next_line
;
; Otherwise, end the segment.
;
finish_segment:
noffset = n_elements(table) / 7
if noffset ne noffset0 then begin
boost_array, iseg, noffset
noffset0 = noffset
endif
;
next_line:
;
endwhile
noffset = n_elements(table) / 7
if noffset ne noffset0 then boost_array, iseg, noffset
free_lun, inlun
;
; Overplot the points in the left and right windows
;
for j=0,n_elements(iseg)-1 do begin
if j eq 0 then i0=0 else i0 = iseg[j-1]
i1 = iseg[j] - 1
setwindow, win_left
plots, table[5,i0:i1], table[6,i0:i1], psym=0, color=color
setwindow, win_right
plots, table[3,i0:i1], table[4,i0:i1], psym=0, color=color
endfor
;
end
;
;==============================================================================
;
pro scc_measure3_replot, left=left, right=right, other=other
;
; Procedure to replot the images, e.g. after changing the image range.
;
common scc_measure3
;
; Unless only the right image was selected, then plot the left image. Get the
; current min and max values.
;
if not keyword_set(right) then begin
if not keyword_set(other) then begin
setwindow, win_left
widget_control, lmin_wid, get_value=lmin1
widget_control, lmax_wid, get_value=lmax1
if lmax1 gt lmin1 then begin
lmin = lmin1
lmax = lmax1
endif
widget_control, lmin_wid, set_value=lmin
widget_control, lmax_wid, set_value=lmax
exptv, subimage_left, /data, /nobox, /noexact, origin=origin_left, $
min=lmin, max=lmax, bscaled=image_left
endif
endif
;
; Unless only the left image was selected, then plot the right image. Get the
; current min and max values.
;
if not keyword_set(left) then begin
if not keyword_set(other) then begin
setwindow, win_right
widget_control, rmin_wid, get_value=rmin1
widget_control, rmax_wid, get_value=rmax1
if rmax1 gt rmin1 then begin
rmin = rmin1
rmax = rmax1
endif
widget_control, rmin_wid, set_value=rmin
widget_control, rmax_wid, set_value=rmax
exptv, subimage_right, /data, /nobox, /noexact, origin=origin_right, $
min=rmin, max=rmax, bscaled=image_right
endif
endif
; Copied for other
;
;
if not keyword_set(left) then begin
if not keyword_set(right) then begin
setwindow, win_other
widget_control, omin_wid, get_value=omin1
widget_control, omax_wid, get_value=omax1
if omax1 gt omin1 then begin
omin = omin1
omax = omax1
endif
widget_control, omin_wid, set_value=omin
widget_control, omax_wid, set_value=omax
exptv, subimage_other, /data, /nobox, /noexact, origin=origin_other, $
min=omin, max=omax, bscaled=image_other
endif
endif
;
; Redisplay any overplots.
;
scc_measure3_redisplay
;
end
;
;==============================================================================
;
pro scc_measure3_select_1st, ev, win_1st, win_2nd, win_3rd,$
image_1st, image_2nd, image_3rd, $
wcs_1st, wcs_2nd, wcs_3rd, sc_1st, sc_2nd, sc_3rd, $
heeq_1st, pix_param_2nd, pix_param_3rd, pix_1st
;
; Procedure to select the first point.
;
common scc_measure3
;
; Start by converting from device pixels to pixels within the image. Overplot
; the selected point.
;
scc_measure3_redisplay, window=win_1st
setwindow, win_1st
pix = convert_coord(ev.x, ev.y, /device, /to_data)
pix_1st = pix[0:1]
if roi[9] then begin ;ROI active
if (pix_1st[0] lt roi[0] or pix_1st[0] gt roi[1] or $
pix_1st[1] lt roi[2] or pix_1st[1] gt roi[3]) then begin
; out of our region of interest, force an early routine
roi[8]=0 ; failed
return
end else begin
roi[8]=1 ; success
endelse
endif
if (sc_1st eq 'ahead') then pixpair[0:1]=pix_1st
if (sc_1st eq 'behind') then pixpair[2:3]=pix_1st
if (sc_1st eq 'eother') then pixpair[4:5]=pix_1st
if (pixforce ne 0) then pixforce=2; mark as 'needs to be saved'
plots, pix_1st[0], pix_1st[1], psym=1, symsize=3, color=color
;
; Convert from image pixel into helioprojective-cartesian coordinates, in
; radians.
;
coord = wcs_get_coord(wcs_1st, pix_1st)
conv = !dpi / 180.d0
case wcs_1st.cunit[0] of
'arcmin': conv = conv / 60.d0
'arcsec': conv = conv / 3600.d0
'mas': conv = conv / 3600.d3
'rad': conv = 1.d0
else: conv = conv
endcase
coord = coord * conv
;
; Calculate the equivalent heliocentric coordinates for distances D within
; +/- maxz of Dsun.
;
dsun = wcs_1st.position.dsun_obs
d = dsun + maxz * [-1,1]
cosy = cos(coord[1])
x = d * cosy * sin(coord[0])
y = d * sin(coord[1])
z = dsun - d * cosy * cos(coord[0])
;
; Determine the spacecraft parameter to pass to convert_sunspice_coord.
;
spacecraft = sc_1st
test = execute('header = h_'+spacecraft)
obsrvtry = 'Earth'
if datatype(header) eq 'STC' then begin
if tag_exist(header, 'OBSRVTRY') then begin
obsrvtry = header.obsrvtry
end else if tag_exist(header, 'TELESCOP') then begin
obsrvtry = header.telescop
endif
end else begin
temp = fxpar(header, 'OBSRVTRY', count=count)
if count gt 0 then obsrvtry = temp else begin
temp = fxpar(header, 'TELESCOP', count=count)
if count gt 0 then obsrvtry = temp
endelse
endelse
spacecraft = parse_sunspice_name(obsrvtry, /earth_default)
if wcs_1st.position.soho and (not wcs_1st.position.pos_assumed) then $
spacecraft = 'SOHO'
IF debugon THEN print,'Selected ',coord,' ',wcs_1st.cunit[0]
if debugon THEN help,dsun,spacecraft
;
; Convert to HEEQ coordinates, with rearranging into HGRTN format as an
; intermediate state.
;
coord = transpose([[z],[x],[y]])
convert_sunspice_coord, wcs_1st.time.observ_date, coord, 'HGRTN', 'HEEQ', $
spacecraft=spacecraft
heeq_1st = coord
;
; Determine the spacecraft parameter to pass to convert_sunspice_coord.
;
spacecraft = sc_2nd
test = execute('header = h_'+spacecraft)
obsrvtry = 'Earth'
if datatype(header) eq 'STC' then begin
if tag_exist(header, 'OBSRVTRY') then begin
obsrvtry = header.obsrvtry
end else if tag_exist(header, 'TELESCOP') then begin
obsrvtry = header.telescop
endif
end else begin
temp = fxpar(header, 'OBSRVTRY', count=count)
if count gt 0 then obsrvtry = temp else begin
temp = fxpar(header, 'TELESCOP', count=count)
if count gt 0 then obsrvtry = temp
endelse
endelse
spacecraft = parse_sunspice_name(obsrvtry, /earth_default)
if wcs_2nd.position.soho and (not wcs_2nd.position.pos_assumed) then $
spacecraft = 'SOHO'
;
; Switch to the other window, and convert from HEEQ to heliocentric-cartesian
; x,y,z values for the other perspective.
;
scc_measure3_redisplay, window=win_2nd
coord2=coord
convert_sunspice_coord, wcs_2nd.time.observ_date, coord2, 'HEEQ', 'HGRTN', $
spacecraft=spacecraft
x = reform(coord2[1,*])
y = reform(coord2[2,*])
z = reform(coord2[0,*])
;
; Convert from heliocentric-cartesian to helioprojective-cartesian. Put into
; the target units.
;
dsun = wcs_2nd.position.dsun_obs
d = sqrt(x^2 + y^2 + (dsun-z)^2)
coord2 = transpose([[atan(x,dsun-z)], [asin(y/d)]])
conv = 180.d0 / !dpi
case wcs_2nd.cunit[0] of
'arcmin': conv = conv * 60.d0
'arcsec': conv = conv * 3600.d0
'mas': conv = conv * 3600.d3
'rad': conv = 1.d0
else: conv = conv
endcase
coord2 = coord2 * conv
IF debugon THEN print,'Translating to ',coord2,' ',wcs_2nd.cunit[0]
IF debugon THEN help,dsun,spacecraft
;
; Convert from helioprojective-cartesian to image pixel coordinates, and
; overplot the points.
;
pix_param_2nd = wcs_get_pixel(wcs_2nd, coord2)
plots, pix_param_2nd[0,*], pix_param_2nd[1,*], color=color
;
; Extrapolate over the full X-range of the image.
;
param_2nd = poly_fit(pix_param_2nd[0,*], pix_param_2nd[1,*], 1)
x = [0,wcs_2nd.naxis[1]-1]
plots, x, poly(x, param_2nd), line=1, color=color
;
; Determine the spacecraft parameter to pass to convert_sunspice_coord.
;
spacecraft = sc_3rd
test = execute('header = h_'+spacecraft)
obsrvtry = 'Earth'
if datatype(header) eq 'STC' then begin
if tag_exist(header, 'OBSRVTRY') then begin
obsrvtry = header.obsrvtry
end else if tag_exist(header, 'TELESCOP') then begin
obsrvtry = header.telescop
endif
end else begin
temp = fxpar(header, 'OBSRVTRY', count=count)
if count gt 0 then obsrvtry = temp else begin
temp = fxpar(header, 'TELESCOP', count=count)
if count gt 0 then obsrvtry = temp
endelse
endelse
spacecraft = parse_sunspice_name(obsrvtry, /earth_default)
if wcs_3rd.position.soho and (not wcs_3rd.position.pos_assumed) then $
spacecraft = 'SOHO'
;
; Switch to the other window, and convert from HEEQ to heliocentric-cartesian
; x,y,z values for the other perspective.
;
scc_measure3_redisplay, window=win_3rd
coord3=coord
convert_sunspice_coord, wcs_3rd.time.observ_date, coord3, 'HEEQ', 'HGRTN', $
spacecraft=spacecraft
x = reform(coord3[1,*])
y = reform(coord3[2,*])
z = reform(coord3[0,*])
;
; Convert from heliocentric-cartesian to helioprojective-cartesian. Put into
; the target units.
;
dsun = wcs_3rd.position.dsun_obs
d = sqrt(x^2 + y^2 + (dsun-z)^2)
coord3 = transpose([[atan(x,dsun-z)], [asin(y/d)]])
conv = 180.d0 / !dpi
case wcs_3rd.cunit[0] of
'arcmin': conv = conv * 60.d0
'arcsec': conv = conv * 3600.d0
'mas': conv = conv * 3600.d3
'rad': conv = 1.d0
else: conv = conv
endcase
coord3 = coord3 * conv
IF debugon THEN print,'Translating to ',coord3,' ',wcs_3rd.cunit[0]
IF debugon THEN help,dsun,spacecraft
;
; Convert from helioprojective-cartesian to image pixel coordinates, and
; overplot the points.
;
pix_param_3rd = wcs_get_pixel(wcs_3rd, coord3)
plots, pix_param_3rd[0,*], pix_param_3rd[1,*], color=color
;
; Extrapolate over the full X-range of the image.
;
param_3rd = poly_fit(pix_param_3rd[0,*], pix_param_3rd[1,*], 1)
x3 = [0,wcs_3rd.naxis[1]-1]
plots, x3, poly(x3, param_3rd), line=1, color=color
;
; Deactivate the zoom button, and clear out the text widgets.
;
widget_control, zoom_wid, sensitive=0
widget_control, out_wid, sensitive=0
widget_control, lon_wid, set_value=''
widget_control, lat_wid, set_value=''
widget_control, rsun_wid, set_value=''
;
end
;
;==============================================================================
;
pro scc_measure3_select_2nd, ev, win_2nd, win_3rd, image_2nd, image_3rd,$
wcs_2nd, wcs_3rd, sc_2nd, sc_3rd, $
heeq_1st, heeq_2nd, pix_param_2nd, pix_param_3rd, $
pix_2nd, pix_3rd
;
; Procedure to select the 2nd point.
;
common scc_measure3
;
; Start by converting from device pixels to pixels within the image.
;
scc_measure3_redisplay, window=win_2nd
pix = convert_coord(ev.x, ev.y, /device, /to_data)
pix = pix[0:1]
; Find the intersection between the linear fit from the other window
; (PARAM_1ST) and the orthogonal line represented by the current selection.
; Overplot the selection.
;
param_2nd = poly_fit(pix_param_2nd[0,*], pix_param_2nd[1,*], 1)
a0 = param_2nd[0]
a1 = param_2nd[1]
x = (pix[0] + a1*(pix[1]-a0)) / (1 + a1^2)
pix_2nd = [x, a0 + a1*x]
;
; Now check our region of interest again
;
if roi[9] then begin ;ROI active
if (pix_2nd[0] lt roi[4] or pix_2nd[0] gt roi[5] or $
pix_2nd[1] lt roi[6] or pix_2nd[1] gt roi[7]) then begin
; out of our region of interest, force an early routine
roi[8]=0 ; failed
return
end else begin
roi[8]=1 ; success
end
endif
if (sc_2nd eq 'ahead') then pixpair[0:1]=pix_2nd
if (sc_2nd eq 'behind') then pixpair[2:3]=pix_2nd
if (sc_2nd eq 'eother') then pixpair[4:5]=pix_2nd
if (pixforce ne 0) then pixforce=2; mark as 'needs to be saved'
plots, pix_2nd[0], pix_2nd[1], psym=1, symsize=3, color=color
;
; Convert from image pixel into helioprojective-cartesian coordinates, in
; radians. Store the angles in the common block
;
coord = wcs_get_coord(wcs_2nd, pix_2nd)
conv = !dpi / 180.d0
case wcs_2nd.cunit[0] of
'arcmin': conv = conv / 60.d0
'arcsec': conv = conv / 3600.d0
'mas': conv = conv / 3600.d3
'rad': conv = 1.d0
else: conv = conv
endcase
coord = coord * conv
;
; Calculate the equivalent heliocentric coordinates for distances D within
; +/- maxz of Dsun.
;
dsun = wcs_2nd.position.dsun_obs
d = dsun + maxz * [-1,1]
cosy = cos(coord[1])
x = d * cosy * sin(coord[0])
y = d * sin(coord[1])
z = dsun - d * cosy * cos(coord[0])
;
; Determine the spacecraft parameter to pass to convert_sunspice_coord.
;
spacecraft = sc_2nd
test = execute('header = h_'+spacecraft)
obsrvtry = 'Earth'
if datatype(header) eq 'STC' then begin
if tag_exist(header, 'OBSRVTRY') then begin
obsrvtry = header.obsrvtry
end else if tag_exist(header, 'TELESCOP') then begin
obsrvtry = header.telescop
endif
end else begin
temp = fxpar(header, 'OBSRVTRY', count=count)
if count gt 0 then obsrvtry = temp else begin
temp = fxpar(header, 'TELESCOP', count=count)
if count gt 0 then obsrvtry = temp
endelse
endelse
spacecraft = parse_sunspice_name(obsrvtry, /earth_default)
if wcs_2nd.position.soho and (not wcs_2nd.position.pos_assumed) then $
spacecraft = 'SOHO'
;
; Convert to HEEQ coordinates, with rearranging into HGRTN format as an
; intermediate state. Store the coordinates in the common block
;
coord = transpose([[z],[x],[y]])
convert_sunspice_coord, wcs_2nd.time.observ_date, coord, 'HGRTN', $
'HEEQ', spacecraft=spacecraft
heeq_2nd = coord
;
spacecraft = sc_3rd
test = execute('header = h_'+spacecraft)
obsrvtry = 'Earth'
if datatype(header) eq 'STC' then begin
if tag_exist(header, 'OBSRVTRY') then begin
obsrvtry = header.obsrvtry
end else if tag_exist(header, 'TELESCOP') then begin
obsrvtry = header.telescop
endif
end else begin
temp = fxpar(header, 'OBSRVTRY', count=count)
if count gt 0 then obsrvtry = temp else begin
temp = fxpar(header, 'TELESCOP', count=count)
if count gt 0 then obsrvtry = temp
endelse
endelse
spacecraft = parse_sunspice_name(obsrvtry, /earth_default)
if wcs_3rd.position.soho and (not wcs_3rd.position.pos_assumed) then $
spacecraft = 'SOHO'
;
; Switch to the other window, and convert from HEEQ to heliocentric-cartesian
; x,y,z values for the other perspective.
; Plot the other line
scc_measure3_redisplay, window=win_3rd
coord3=coord
convert_sunspice_coord, wcs_3rd.time.observ_date, coord3, 'HEEQ', 'HGRTN', $
spacecraft=spacecraft
x = reform(coord3[1,*])
y = reform(coord3[2,*])
z = reform(coord3[0,*])
;
; Convert from heliocentric-cartesian to helioprojective-cartesian. Put into
; the target units.
;
dsun = wcs_3rd.position.dsun_obs
d = sqrt(x^2 + y^2 + (dsun-z)^2)
coord3 = transpose([[atan(x,dsun-z)], [asin(y/d)]])
conv = 180.d0 / !dpi
case wcs_3rd.cunit[0] of
'arcmin': conv = conv * 60.d0
'arcsec': conv = conv * 3600.d0
'mas': conv = conv * 3600.d3
'rad': conv = 1.d0
else: conv = conv
endcase
coord3 = coord3 * conv
IF debugon THEN print,'Translating to ',coord3,' ',wcs_3rd.cunit[0]
IF debugon THEN help,dsun,spacecraft
;
; Convert from helioprojective-cartesian to image pixel coordinates, and
; overplot the points.
;
plots, pix_param_3rd[0,*], pix_param_3rd[1,*], color=color
pix_param4 = wcs_get_pixel(wcs_3rd, coord3)
plots, pix_param4[0,*], pix_param4[1,*], color=color1
;
; Extrapolate over the full X-range of the image.
;
param_4 = poly_fit(pix_param4[0,*], pix_param4[1,*], 1)
x3 = [0,wcs_3rd.naxis[1]-1]
plots, x3, poly(x3, param_4), line=1, color=color1
;
; Based on the HEEQ coordinates from the left and right images, find the
; intersection on the equatorial (x-y) plane.
;
p1st = poly_fit(heeq_1st[0,*], heeq_1st[1,*], 1)
p2nd = poly_fit(heeq_2nd[0,*], heeq_2nd[1,*], 1)
x = (p1st[0] - p2nd[0]) / (p2nd[1] - p1st[1])
y = (poly(x,p1st) + poly(x,p2nd)) / 2
;
; Using the same point, find the Z position.
;
p1st = poly_fit(heeq_1st[0,*], heeq_1st[2,*], 1)
p2nd = poly_fit(heeq_2nd[0,*], heeq_2nd[2,*], 1)
z = (poly(x,p1st) + poly(x,p2nd)) / 2
;
; Store 3rd image coordinates of the point
coordhq = transpose([[x],[y],[z]])
convert_sunspice_coord, wcs_3rd.time.observ_date, coordhq, 'HEEQ', 'HGRTN', $
spacecraft=spacecraft
x3 = reform(coordhq[1,*])
y3 = reform(coordhq[2,*])
z3 = reform(coordhq[0,*])
;
; Convert from heliocentric-cartesian to helioprojective-cartesian. Put into
; the target units.
;
dsun = wcs_3rd.position.dsun_obs
d = sqrt(x3^2 + y3^2 + (dsun-z3)^2)
coordhq = transpose([[atan(x3,dsun-z3)], [asin(y3/d)]])
conv = 180.d0 / !dpi
case wcs_3rd.cunit[0] of
'arcmin': conv = conv * 60.d0
'arcsec': conv = conv * 3600.d0
'mas': conv = conv * 3600.d3
'rad': conv = 1.d0
else: conv = conv
endcase
coordhq= coordhq * conv
;
; Convert from helioprojective-cartesian to image pixel coordinates, and
; overplot the points.
;
pix_3rd = wcs_get_pixel(wcs_3rd, coordhq)
; Now check our region of interest again
;
if roi[9] then begin ;ROI active
if (pix_3rd[0] lt roi[10] or pix_3rd[0] gt roi[11] or $
pix_3rd[1] lt roi[12] or pix_3rd[1] gt roi[13]) then begin
; out of our region of interest, force an early routine
roi[8]=0 ; failed
return
end else begin
roi[8]=1 ; success
end
endif
if (sc_3rd eq 'ahead') then pixpair[0:1]=pix_3rd
if (sc_3rd eq 'behind') then pixpair[2:3]=pix_3rd
if (sc_3rd eq 'eother') then pixpair[4:5]=pix_3rd
plots, pix_3rd[0], pix_3rd[1], psym=2, symsize=3, color=color2
;
; Populate the widgets.
;
rad = sqrt(x^2 + y^2 + z^2)
lon = atan(y, x) * 180 / !dpi
widget_control, lon_wid, set_value=ntrim(float(lon))
lat = asin(z / rad) * 180 / !dpi
widget_control, lat_wid, set_value=ntrim(float(lat))
rad = rad / 6.95508e8
widget_control, rsun_wid, set_value=ntrim(float(rad))
;
; Activate the zoom and store buttons.
;
widget_control, zoom_wid, /sensitive
if outlun gt 0 then widget_control, out_wid, /sensitive
;
end
;
;==============================================================================
;
pro scc_measure3_cleanup, id ;without modifications from 2img version
;
; Cleanup procedure, to make sure that the output file is properly closed.
;
common scc_measure3
if outlun gt 0 then free_lun, outlun
end
;
;==============================================================================
;
pro store_scc_measure3
common scc_measure3
widget_control, lon_wid, get_value=lon
widget_control, lat_wid, get_value=lat
widget_control, rsun_wid, get_value=rad
printf, outlun, lon, lat, rad, pixpair, format='(9(F11.5,1x))'
; printf, outlun, lon, lat, rad, format='(3F5.5)'
flush, outlun
widget_control, out_wid, sensitive=0
widget_control, clear_wid, /sensitive
widget_control, oplot_wid, /sensitive
n_stored = n_stored + 1
if n_stored eq 1 then begin
store_left = pix_left
store_right = pix_right
store_other = pix_other
end else begin
boost_array, store_left, pix_left
boost_array, store_right, pix_right
boost_array, store_other, pix_other
endelse
end
;==============================================================================
;
pro scc_measure3_event, ev
;
; Widget event handler.
;
common scc_measure3
;
widget_control, ev.id, get_uvalue=uvalue
case uvalue of
'EXIT': begin
if (pixforce eq 2) then store_scc_measure3
widget_control, /destroy, ev.top
end
;
'XLOADCT': xloadct, group=ev.top, updatecallback="scc_measure3_redisplay"
;
; The plot color was changed.
;
'COLOR': begin
widget_control, ev.id, get_value=color
color = 0 > color < (!d.table_size - 1)
widget_control, ev.id, set_value=color
scc_measure3_replot
end
;
; The left window was selected. Start by converting from device pixels to
; pixels within the image. Overplot the selected point.
;
'LEFT': if ev.press gt 0 then begin
if win_last eq 'LEFT' then in_progress = 0
if in_progress then begin
if win_last eq 'RIGHT' then begin
scc_measure3_select_2nd, ev, win_left, win_other, $
image_left, image_other, wcs_behind, wcs_eother, $
'behind','eother', heeq_right, heeq_left, $
pix_paraml, pix_paramo, pix_left, pix_other
if roi[8] or (roi[9] eq 0) then in_progress = 0
end else begin
scc_measure3_select_2nd, ev, win_left, win_right, $
image_left, image_right, wcs_behind, wcs_ahead, $
'behind','ahead', heeq_other, heeq_left, $
pix_paraml, pix_paramr, pix_left, pix_right
if roi[8] or (roi[9] eq 0) then in_progress = 0
endelse
end else begin
scc_measure3_select_1st, ev, win_left, win_right, win_other, $
image_left, image_right, image_other, wcs_behind, $
wcs_ahead, wcs_eother, 'behind', 'ahead', 'eother', $
heeq_left, pix_paramr, pix_paramo, pix_left
if roi[8] or (roi[9] eq 0) then in_progress = 1
endelse
if roi[8] or (roi[9] eq 0) then win_last = 'LEFT'
endif
;
; The right window was selected. Start by converting from device pixels to
; pixels within the image.
;
'RIGHT': if ev.press gt 0 then begin
if win_last eq 'RIGHT' then in_progress = 0
if in_progress then begin
if win_last eq 'LEFT' then begin
scc_measure3_select_2nd, ev, win_right, win_other, $
image_right, image_other, wcs_ahead, wcs_eother, $
'ahead','eother', heeq_left, heeq_right, $
pix_paramr, pix_paramo, pix_right, pix_other
if roi[8] or (roi[9] eq 0) then in_progress = 0
end else begin
scc_measure3_select_2nd, ev, win_right, win_left, $
image_right, image_left, wcs_ahead, wcs_behind, $
'ahead','behind', heeq_other, heeq_right, $
pix_paramr, pix_paraml, pix_right, pix_left
if roi[8] or (roi[9] eq 0) then in_progress = 0
endelse
end else begin
scc_measure3_select_1st, ev, win_right, win_left, win_other, $
image_right, image_left, image_other, wcs_ahead,$
wcs_behind, wcs_eother, 'ahead', 'behind', 'eother', $
heeq_right, pix_paraml, pix_paramo, pix_right
if roi[8] or (roi[9] eq 0) then in_progress = 1
endelse
if roi[8] or (roi[9] eq 0) then win_last = 'RIGHT'
endif
;
; The right window was selected. Start by converting from device pixels to
; pixels within the image.
;
'OTHER': if ev.press gt 0 then begin
if win_last eq 'OTHER' then in_progress = 0
if in_progress then begin
if win_last eq 'RIGHT' then begin
scc_measure3_select_2nd, ev, win_other, win_left, $
image_other, image_left, wcs_eother, wcs_behind, $
'eother','behind', heeq_right, heeq_other, $
pix_paramo, pix_paraml, pix_other, pix_left
if roi[8] or (roi[9] eq 0) then in_progress = 0
end else begin
scc_measure3_select_2nd, ev, win_other, win_right, $
image_other, image_right, wcs_eother, wcs_ahead, $
'eother','ahead', heeq_left, heeq_other, $
pix_paramo, pix_paramr, pix_other, pix_right
if roi[8] or (roi[9] eq 0) then in_progress = 0
endelse
end else begin
scc_measure3_select_1st, ev, win_other, win_right, win_left, $
image_other, image_right, image_left, wcs_eother,$
wcs_ahead, wcs_behind, 'eother', 'ahead', 'behind', $
heeq_other, pix_paramr, pix_paraml, pix_other
if roi[8] or (roi[9] eq 0) then in_progress = 1
endelse
if roi[8] or (roi[9] eq 0) then win_last = 'OTHER'
endif
;
; Modify the image ranges.
;
'LMIN': scc_measure3_replot, /left
'LMAX': scc_measure3_replot, /left
'RMIN': scc_measure3_replot, /right
'RMAX': scc_measure3_replot, /right
'OMIN': scc_measure3_replot, /other
'OMAX': scc_measure3_replot, /other
;
; Store any selected points.
;
'STORE': begin
store_scc_measure3
if (pixforce ne 0) then pixforce=1; mark as 'saved'
end
;
; Clear the previously stored points from memory.
;
'CLEAR': begin
n_stored = 0
printf, outlun
flush, outlun
widget_control, clear_wid, sensitive=0
scc_measure3_redisplay
end
;
; Overplot previously stored points from the output file.
;
'OPLOT': scc_measure3_oplot_from_file
;
; Zoom in on the previously selected points.
;
'ZOOM': begin