% clear workspace
clear all

% close all figures
close all

% -----------------------------------------------------------------------------

cd C:\SOFIA_Octave
pause(1)

dir
pause(1)

pkg list fits
pause(1)

pkg load fits
pause(1)

pkg describe fits -verbose
pause(1)

read_fits_image("F0517_EX_SPE_0601176_EXEELONEXEECHL_MRD_0023-0029.fits")
pause(1)

[DATA, HEADER] = read_fits_image("F0517_EX_SPE_0601176_EXEELONEXEECHL_MRD_0023-0029.fits")
pause(1)

whos
pause(1)

HEADER
pause(1)

plot( DATA(:,2) )
pause(1)

plot( DATA(:, 1), DATA(:,2) )
pause(1)

plot( 10000 ./ DATA(:, 1), DATA(:,2) )
pause(1)

lambda = 10000 ./ DATA(:, 1)
pause(1)

intensity_instrument = DATA(:,2)
pause(1)

plot(lambda, intensity_instrument)
pause(1)

plot(lambda, intensity_instrument, '-k')
pause(1)

grid on
pause(1)

atmospheric_transmission = DATA(:,4)
pause(1)

plot(lambda, atmospheric_transmission, '-r')
grid on
pause(1)

plot( lambda, intensity_instrument, '-k', lambda, atmospheric_transmission, '-r' )
grid on
pause(1)

max(intensity_instrument)
pause(1)

plot( lambda, intensity_instrument, '-k', ...
lambda, atmospheric_transmission * max(intensity_instrument), '-r' )
grid on
pause(1)

intensity_instrument
pause(1)

intensity_instrument < 0
pause(1)

idx = find( intensity_instrument < 0 )
pause(1)

size(idx)
pause(1)

plot( lambda, intensity_instrument, '-k', lambda(idx), intensity_instrument(idx), '+g' )
grid on
pause(1)

intensity_instrument(idx) = 0
pause(1)

plot( lambda, intensity_instrument, '-k', ...
lambda, atmospheric_transmission * max(intensity_instrument), '-r' )
grid on
pause(1)

xlabel("Wellenlaenge [Mikrometer]")
ylabel("Energiefluss [erg / (s cm2 sr (cm-1)]")
title("SOFIA-EXES-Spektrum des Protosterns AFGL 2136")
pause(1)

% -----------------------------------------------------------------------------
