hrt_pipe.py 35.2 KB
Newer Older
1
from posix import listdir
jonas's avatar
jonas committed
2
3
4
5
import numpy as np 
import os.path
from astropy.io import fits
from scipy.ndimage import gaussian_filter
6
import time
jonas's avatar
jonas committed
7
from operator import itemgetter
jonas's avatar
jonas committed
8
import json
jonas's avatar
jonas committed
9
import matplotlib.pyplot as plt
jonas's avatar
jonas committed
10

Jonas Sinjan's avatar
Jonas Sinjan committed
11
from utils import *
jonas's avatar
jonas committed
12
13


jonas's avatar
jonas committed
14
def get_data(path, scaling = True, bit_convert_scale = True):
jonas's avatar
jonas committed
15
    """load science data from path"""
jonas's avatar
jonas committed
16
17
    #try:
    data, header = load_fits(path)
18

jonas's avatar
jonas committed
19
20
    if bit_convert_scale: #conversion from 24.8bit to 32bit
        data /=  256.
jonas's avatar
jonas committed
21

jonas's avatar
jonas committed
22
23
24
    if scaling:
        
        accu = header['ACCACCUM']*header['ACCROWIT']*header['ACCCOLIT'] #getting the number of accu from header
jonas's avatar
jonas committed
25

jonas's avatar
jonas committed
26
        data /= accu
jonas's avatar
jonas committed
27

jonas's avatar
jonas committed
28
29
30
        printc(f"Dividing by number of accumulations: {accu}",color=bcolors.OKGREEN)
    
    return data, header
jonas's avatar
jonas committed
31

jonas's avatar
jonas committed
32
33
    #except Exception:
    #    printc("ERROR, Unable to open fits file: {}",path,color=bcolors.FAIL)
jonas's avatar
jonas committed
34

Jonas Sinjan's avatar
Jonas Sinjan committed
35
   
jonas's avatar
jonas committed
36

Jonas Sinjan's avatar
Jonas Sinjan committed
37
def demod_hrt(data,pmp_temp):
jonas's avatar
jonas committed
38
    '''
Jonas Sinjan's avatar
Jonas Sinjan committed
39
    Use constant demodulation matrices to demodulate data
jonas's avatar
jonas committed
40
    '''
Jonas Sinjan's avatar
Jonas Sinjan committed
41
42
43
44
45
46
 
    if pmp_temp == '50':
        demod_data = np.array([[ 0.28037298,  0.18741922,  0.25307596,  0.28119895],
                     [ 0.40408596,  0.10412157, -0.7225681,   0.20825675],
                     [-0.19126636, -0.5348939,   0.08181918,  0.64422774],
                     [-0.56897295,  0.58620095, -0.2579202,   0.2414017 ]])
Jonas Sinjan's avatar
Jonas Sinjan committed
47
        
Jonas Sinjan's avatar
Jonas Sinjan committed
48
49
50
51
52
    elif pmp_temp == '40':
        demod_data = np.array([[ 0.26450154,  0.2839626,   0.12642948,  0.3216773 ],
                     [ 0.59873885,  0.11278069, -0.74991184,  0.03091451],
                     [ 0.10833212, -0.5317737,  -0.1677862,   0.5923593 ],
                     [-0.46916953,  0.47738808, -0.43824592,  0.42579797]])
jonas's avatar
jonas committed
53
    
jonas's avatar
jonas committed
54
    else:
Jonas Sinjan's avatar
Jonas Sinjan committed
55
        printc("Demodulation Matrix for PMP TEMP of {pmp_temp} deg is not available", color = bcolors.FAIL)
jonas's avatar
jonas committed
56

Jonas Sinjan's avatar
Jonas Sinjan committed
57
    printc(f'Using a constant demodulation matrix for a PMP TEMP of {pmp_temp} deg',color = bcolors.OKGREEN)
jonas's avatar
jonas committed
58
    
Jonas Sinjan's avatar
Jonas Sinjan committed
59
60
61
    demod_data = demod_data.reshape((4,4))
    shape = data.shape
    demod = np.tile(demod_data, (shape[0],shape[1],1,1))
jonas's avatar
jonas committed
62
63
64
65

    if data.ndim == 5:
        #if data array has more than one scan
        data = np.moveaxis(data,-1,0) #moving number of scans to first dimension
jonas's avatar
jonas committed
66

jonas's avatar
jonas committed
67
68
        data = np.matmul(demod,data)
        data = np.moveaxis(data,0,-1) #move scans back to the end
jonas's avatar
jonas committed
69
70
71
    
    elif data.ndim == 4:
        #for if data has just one scan
jonas's avatar
jonas committed
72
        data = np.matmul(demod,data)
jonas's avatar
jonas committed
73
    
jonas's avatar
jonas committed
74
    return data, demod
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

def check_filenames(data_f):
    #checking if the science scans have the same DID - this would cause an issue for naming the output demod files

    scan_name_list = [str(scan.split('.fits')[0][-10:]) for scan in data_f]

    seen = set()
    uniq_scan_DIDs = [x for x in scan_name_list if x in seen or seen.add(x)] #creates list of unique DIDs from the list

    #print(uniq_scan_DIDs)
    #print(scan_name_list)
    if uniq_scan_DIDs == []:
        print("The scans' DIDs are all unique")

    else:

        for x in uniq_scan_DIDs:
            number = scan_name_list.count(x)
            if number > 1: #if more than one
                print(f"The DID: {x} is repeated {number} times")
                i = 1
                for index, name in enumerate(scan_name_list):
                    if name == x:
                        scan_name_list[index] = name + f"_{i}" #add _1, _2, etc to the file name, so that when written to output file not overwriting
                        i += 1

        print("The New DID list is: ", scan_name_list)

    return scan_name_list
jonas's avatar
jonas committed
104
105
106
    


107
def phihrt_pipe(data_f, dark_f = '', flat_f = '', L1_input = True, L1_8_generate = False, scale_data = True, accum_scaling = True, 
jonas's avatar
jonas committed
108
                bit_conversion = True, norm_f = True, clean_f = False, sigma = 59, clean_mode = "V", flat_states = 24, prefilter_f = None,flat_c = True, 
109
                dark_c = True, fs_c = True, demod = True, norm_stokes = True, out_dir = './',  out_demod_file = False,  out_demod_filename = None,
jonas's avatar
jonas committed
110
                ItoQUV = False, ctalk_params = None, rte = False, out_rte_filename = None, p_milos = True, config_file = True):
jonas's avatar
jonas committed
111
112
113

    '''
    PHI-HRT data reduction pipeline
jonas's avatar
jonas committed
114
    1. read in science data (+ OPTION: scaling) open path option + open for several scans at once
jonas's avatar
jonas committed
115
116
117
118
    2. read in flat field (+scaling)- just accepts one flat field fits file
    3. read in dark field (+scaling)
    4. apply dark field
    5. option to clean flat field with unsharp masking
jonas's avatar
jonas committed
119
120
    6. normalise flat field
    7. apply flat field
jonas's avatar
jonas committed
121
    8. apply prefilter
jonas's avatar
jonas committed
122
123
124
    9. read in field stop
    10. apply field stop
    11. demodulate with const demod matrix
jonas's avatar
jonas committed
125
        a) option to output demod to fits file
jonas's avatar
jonas committed
126
127
    12. normalise to quiet sun
    13. calibration
jonas's avatar
jonas committed
128
        a) ghost correction - not implemented yet
jonas.sinjan's avatar
jonas.sinjan committed
129
        b) cross talk correction
jonas's avatar
jonas committed
130
    14. rte inversion with pmilos/cmilos (CE, RTE or CE+RTE)
jonas's avatar
jonas committed
131
        a) output rte data products to fits file
jonas's avatar
jonas committed
132
133
134
135
136

    Parameters
    ----------
        Input:
    data_f : list or string
jonas's avatar
jonas committed
137
        list containing paths to fits files of the raw HRT data OR string of path to one file  - must have last 10 characters before .fits as the DID - for naming purposes of output files
138
    dark_f : string, DEFAULT ''
jonas's avatar
jonas committed
139
        Fits file of a dark file (ONLY ONE FILE)
140
    flat_f : string, DEFAULT ''
jonas's avatar
jonas committed
141
142
143
        Fits file of a HRT flatfield (ONLY ONE FILE)

    ** Options:
144
145
146
147
    L1_input: bool, DEFAULT True
        ovverides scale_data, bit_conversion, and accum_scaling, so that correct scaling for L1 data applied
    L1_8_generate: bool, DEFAULT False
        if True, assumes L1 input, and generates RTE output with the calibration header information
148
149
    scale_data: bool, DEFAULT True
        performs the accumulation scaling + conversion for flat and science (only FALSE for commissioning data)
150
    bit_conversion: bool, DEFAULT True
151
        divides the scan + flat by 256 to convert from 24.8bit to 32bits
jonas's avatar
jonas committed
152
153
154
    norm_f: bool, DEFAULT: True
        to normalise the flat fields before applying
    clean_f: bool, DEFAULT: False
jonas's avatar
jonas committed
155
156
        clean the flat field with unsharp masking
    sigma: int, DEFAULT: 59
jonas's avatar
jonas committed
157
158
159
        sigma of the gaussian convolution used for unsharp masking if clean_f == True
    clean_mode: str, DEFAULT: "V"
        The polarisation states of the flat field to be unsharp masked, options are "V", "UV" and "QUV"
jonas's avatar
jonas committed
160
161
    flat_states: int, DEFAULT: 24
        Number of flat fields to be applied, options are 4 (one for each pol state), 6 (one for each wavelength), 24 (one for each image)
jonas's avatar
jonas committed
162
163
    prefilter_f: str, DEFAULT None
        file path location to prefilter fits file, apply prefilter correction
jonas's avatar
jonas committed
164
165
166
167
    flat_c: bool, DEFAULT: True
        apply flat field correction
    dark_c: bool, DEFAULT: True
        apply dark field correction
jonas's avatar
jonas committed
168
    fs_c: bool, DEFAULT True
169
        apply HRT field stop
jonas's avatar
jonas committed
170
171
172
173
174
175
176
177
    demod: bool, DEFAULT: True
        apply demodulate to the stokes
    norm_stokes: bool, DEFAULT: True
        normalise the stokes vector to the quiet sun (I_continuum)
    out_dir : string, DEFUALT: './'
        directory for the output files
    out_demod_file: bool, DEFAULT: False
        output file with the stokes vectors to fits file
178
179
    out_demod_filename: str, DEFAULT = None
        if None, takes last 10 characters of input scan filename (assumes its a DID), change if want other name
jonas's avatar
jonas committed
180
181
    ItoQUV: bool, DEFAULT: False 
        apply I -> Q,U,V correction
jonas's avatar
jonas committed
182
183
    ctalk_params: numpy arr, DEFAULT: None 
        cross talk parameters for ItoQUV, (2,3) numpy array required: first axis: Slope, Offset (Normalised to I_c) - second axis:  Q,U,V
jonas's avatar
jonas committed
184
185
    rte: str, DEFAULT: False 
        invert using cmilos, options: 'RTE' for Milne Eddington Inversion, 'CE' for Classical Estimates, 'CE+RTE' for combined
186
187
    out_rte_filename: str, DEFAULT = ''
        if '', takes last 10 characters of input scan filename (assumes its a DID), change if want other name
188
    p_milos: bool, DEFAULT = True
jonas's avatar
jonas committed
189
        if True, will execute the RTE inversion using the parallel version of the CMILOS code on 16 processors
jonas's avatar
jonas committed
190
191
    config_file: bool, DEFAULT = True
        if True, will generate config.txt file that writes the reduction process steps done
jonas's avatar
jonas committed
192
193
194
    
    Returns
    -------
jonas's avatar
jonas committed
195
196
    data: nump array
        stokes vector 
jonas's avatar
jonas committed
197
198
199
200
201
202
203
204
205
206
207

    References
    ----------
    SPGYlib

    '''

    printc('--------------------------------------------------------------',bcolors.OKGREEN)
    printc('PHI HRT data reduction software  ',bcolors.OKGREEN)
    printc('--------------------------------------------------------------',bcolors.OKGREEN)

Jonas Sinjan's avatar
Jonas Sinjan committed
208
    overall_time = time.time()
jonas's avatar
jonas committed
209
    saved_args = locals()
jonas's avatar
jonas committed
210
    saved_args['ctalk_params'] = ctalk_params.tolist()
211
212
213
214
215
216
217
218


    if L1_input:
        print("L1_input param set to True - Assuming L1 science data")
        accum_scaling = True
        bit_conversion = False
        scale_data = True

jonas's avatar
jonas committed
219
220
221
    #-----------------
    # READ DATA
    #-----------------
222

Jonas Sinjan's avatar
Jonas Sinjan committed
223
    print(" ")
Jonas Sinjan's avatar
Jonas Sinjan committed
224
    printc('-->>>>>>> Reading Data              ',color=bcolors.OKGREEN) 
jonas's avatar
jonas committed
225

226
    start_time = time.time()
jonas's avatar
jonas committed
227

228
229
230
    if isinstance(data_f, str):
        data_f = [data_f]

jonas's avatar
jonas committed
231
    if isinstance(data_f, list):
232
        #if the data_f contains several scans
Jonas Sinjan's avatar
Jonas Sinjan committed
233
        printc(f'Input contains {len(data_f)} scan(s)',color=bcolors.OKGREEN)
234
235
236
237
238
239
        
        number_of_scans = len(data_f)

        data_arr = [0]*number_of_scans
        hdr_arr = [0]*number_of_scans

jonas's avatar
jonas committed
240
241
        wve_axis_arr = [0]*number_of_scans
        cpos_arr = [0]*number_of_scans
jonas's avatar
jonas committed
242
243
        voltagesData_arr = [0]*number_of_scans
        tuning_constant_arr = [0]*number_of_scans
jonas's avatar
jonas committed
244

245
        for scan in range(number_of_scans):
246
            data_arr[scan], hdr_arr[scan] = get_data(data_f[scan], scaling = accum_scaling, bit_convert_scale = bit_conversion)
247

248
            wve_axis_arr[scan], voltagesData_arr[scan], tuning_constant_arr[scan], cpos_arr[scan] = fits_get_sampling(data_f[scan],verbose = True)
249

jonas's avatar
jonas committed
250
            if scale_data: #not for commissioning data
251

jonas's avatar
jonas committed
252
                data_arr[scan] *= 81920/127 #conversion factor if 16 bits
jonas's avatar
jonas committed
253

jonas's avatar
jonas committed
254
            if 'IMGDIRX' in hdr_arr[scan] and hdr_arr[scan]['IMGDIRX'] == 'YES':
255
                print(f"This scan has been flipped in the Y axis to conform to orientation standards. \n File: {data_f[scan]}")
jonas's avatar
jonas committed
256

jonas's avatar
jonas committed
257

jonas's avatar
jonas committed
258
        #--------
jonas's avatar
jonas committed
259
        # test if the scans have different sizes
jonas's avatar
jonas committed
260
261
        #--------

262
263
264
        first_shape = data_arr[scan].shape
        result = all(element.shape == first_shape for element in data_arr)
        if (result):
265
            print("All the scan(s) have the same dimension")
jonas's avatar
jonas committed
266

267
268
269
270
271
        else:
            print("The scans have different dimensions! \n Ending process")

            exit()

jonas's avatar
jonas committed
272
273

        #--------
jonas's avatar
jonas committed
274
        # test if the scans have different continuum wavelength_positions
jonas's avatar
jonas committed
275
276
277
278
279
        #--------

        first_cpos = cpos_arr[0]
        result = all(c_position == first_cpos for c_position in cpos_arr)
        if (result):
280
            print("All the scan(s) have the same continuum wavelength position")
jonas's avatar
jonas committed
281
282

        else:
283
284
285
            print("The scans have different continuum_wavelength postitions! Please fix \n Ending Process")

            exit()
jonas's avatar
jonas committed
286

jonas's avatar
jonas committed
287
        #--------
jonas's avatar
jonas committed
288
        # test if the scans have different pmp temperatures
jonas's avatar
jonas committed
289
290
291
292
293
        #--------

        first_pmp_temp = hdr_arr[0]['HPMPTSP1']
        result = all(hdr['HPMPTSP1'] == first_pmp_temp for hdr in hdr_arr)
        if (result):
294
            print(f"All the scan(s) have the same PMP Temperature Set Point: {first_pmp_temp}")
jonas's avatar
jonas committed
295
296
297
298
299
300
301
            pmp_temp = str(first_pmp_temp)

        else:
            print("The scans have different PMP Temperatures! Please fix \n Ending Process")

            exit()

302
303
        data = np.stack(data_arr, axis = -1)
        data = np.moveaxis(data, 0,-2) #so that it is [y,x,24,scans]
jonas's avatar
jonas committed
304

305
        print(f"Data shape is {data.shape}")
jonas's avatar
jonas committed
306

jonas's avatar
jonas committed
307
308
309
310
        #--------
        # test if the scans have same IMGDIRX keyword
        #--------
        if all('IMGDIRX' in hdr for hdr in hdr_arr):
jonas's avatar
jonas committed
311
            header_imgdirx_exists = True
jonas's avatar
jonas committed
312
313
314
            first_imgdirx = hdr_arr[0]['IMGDIRX']
            result = all(hdr['IMGDIRX'] == first_imgdirx for hdr in hdr_arr)
            if (result):
315
                print(f"All the scan(s) have the same IMGDIRX keyword: {first_imgdirx}")
jonas's avatar
jonas committed
316
317
318
319
320
                imgdirx_flipped = str(first_imgdirx)
            else:
                print("The scans have different IMGDIRX keywords! Please fix \n Ending Process")
                exit()
        else:
jonas's avatar
jonas committed
321
            header_imgdirx_exists = False
322
            print("Not all the scan(s) contain the 'IMGDIRX' keyword! Assuming all not flipped - Proceed with caution")
jonas's avatar
jonas committed
323

jonas's avatar
jonas committed
324
    else:
325
326
327
        printc("ERROR, data_f argument is neither a string nor list containing strings: {} \n Ending Process",data_f,color=bcolors.FAIL)
        exit()

jonas's avatar
jonas committed
328
329
    saved_args['science_cpos'] = int(cpos_arr[0])

330
331
332
333
334
335
336
337
338
339
    data_shape = data.shape

    data_size = data_shape[:2]
    
    #converting to [y,x,pol,wv,scans]

    data = data.reshape(data_size[0],data_size[1],6,4,data_shape[-1]) #separate 24 images, into 6 wavelengths, with each 4 pol states
    data = np.moveaxis(data, 2,-2) #need to swap back to work

    #enabling cropped datasets, so that the correct regions of the dark field and flat field are applied
jonas.sinjan's avatar
jonas.sinjan committed
340
    print("Data reshaped to: ", data.shape)
341

Jonas Sinjan's avatar
Jonas Sinjan committed
342
    diff = 2048-data_size[0] #handling 0/2 errors
343
344
345
346
347
348
349
350
    
    if np.abs(diff) > 0:
    
        start_row = int((2048-data_size[0])/2)
        start_col = int((2048-data_size[1])/2)
        
    else:
        start_row, start_col = 0, 0
Jonas Sinjan's avatar
Jonas Sinjan committed
351
352
353
354
        
    printc('--------------------------------------------------------------',bcolors.OKGREEN)
    printc(f"------------ Load science data time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
    printc('--------------------------------------------------------------',bcolors.OKGREEN)
jonas.sinjan's avatar
jonas.sinjan committed
355

jonas's avatar
jonas committed
356
357
358
359
360
    #-----------------
    # READ FLAT FIELDS
    #-----------------

    if flat_c:
Jonas Sinjan's avatar
Jonas Sinjan committed
361
        print(" ")
362
        printc('-->>>>>>> Reading Flats',color=bcolors.OKGREEN)
jonas's avatar
jonas committed
363

364
365
        start_time = time.time()
    
jonas's avatar
jonas committed
366
        #try:
367
        flat, header_flat = get_data(flat_f, scaling = accum_scaling,  bit_convert_scale=bit_conversion)
368

jonas's avatar
jonas committed
369
370
371
        print(f"Flat field shape is {flat.shape}")
        
        if header_flat['BITPIX'] == 16:
372

jonas's avatar
jonas committed
373
            print("Number of bits per pixel is: 16")
374

jonas's avatar
jonas committed
375
            flat *= 614400/128
jonas's avatar
jonas committed
376

jonas's avatar
jonas committed
377
378
379
380
381
382
383
        flat = np.moveaxis(flat, 0,-1) #so that it is [y,x,24]
        flat = flat.reshape(2048,2048,6,4) #separate 24 images, into 6 wavelengths, with each 4 pol states
        flat = np.moveaxis(flat, 2,-1)
        
        print(flat.shape)

        _, _, _, cpos_f = fits_get_sampling(flat_f,verbose = True)
jonas's avatar
jonas committed
384

jonas's avatar
jonas committed
385
        print(f"The continuum position of the flat field is at {cpos_f} index position")
jonas's avatar
jonas committed
386
387

        saved_args['flat_cpos'] = int(cpos_f)
jonas's avatar
jonas committed
388
        
jonas's avatar
jonas committed
389
390
        if cpos_f != cpos_arr[0]:
            print("The flat field continuum position is not the same as the data, trying to correct.")
391

jonas's avatar
jonas committed
392
            if cpos_f == 5 and cpos_arr[0] == 0:
jonas's avatar
jonas committed
393

jonas's avatar
jonas committed
394
                flat = np.roll(flat, 1, axis = -1)
jonas's avatar
jonas committed
395

jonas's avatar
jonas committed
396
397
398
399
400
401
402
403
            elif cpos_f == 0 and cpos_arr[0] == 5:

                flat = np.roll(flat, -1, axis = -1)

            else:
                print("Cannot reconcile the different continuum positions. \n Ending Process.")

                exit()
jonas's avatar
jonas committed
404

jonas's avatar
jonas committed
405
        flat_pmp_temp = str(header_flat['HPMPTSP1'])
jonas's avatar
jonas committed
406

jonas's avatar
jonas committed
407
        print(f"Flat PMP Temperature Set Point: {flat_pmp_temp}")
408
409

        if flat_f[-62:] == 'solo_L0_phi-hrt-flat_0667134081_V202103221851C_0162201100.fits': 
jonas's avatar
jonas committed
410
411
412
413
            print("This flat has a missing line - filling in with neighbouring pixels")
            flat_copy = flat.copy()
            flat[1345, 296:, 1, 1] = flat_copy[1344, 296:, 1, 1]
            flat[1346, :291, 1, 1] = flat_copy[1345, :291, 1, 1]
jonas's avatar
jonas committed
414
415
416
417
            
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------ Load flats time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
418

jonas's avatar
jonas committed
419
420
        # except Exception:
        #     printc("ERROR, Unable to open/process flats file: {}",flat_f,color=bcolors.FAIL)
jonas's avatar
jonas committed
421
422
423


    else:
Jonas Sinjan's avatar
Jonas Sinjan committed
424
        print(" ")
425
426
        printc('-->>>>>>> No flats mode',color=bcolors.WARNING)

jonas's avatar
jonas committed
427
428
429
430

    #-----------------
    # READ AND CORRECT DARK FIELD
    #-----------------
431

jonas's avatar
jonas committed
432
    if dark_c:
Jonas Sinjan's avatar
Jonas Sinjan committed
433
        print(" ")
jonas's avatar
jonas committed
434
        printc('-->>>>>>> Reading Darks                   ',color=bcolors.OKGREEN)
435
436
437

        start_time = time.time()

jonas's avatar
jonas committed
438
        try:
jonas's avatar
jonas committed
439
440
441
442
443
444
445
446
447
448
449
450
451
452
            dark,h = get_data(dark_f)

            dark_shape = dark.shape

            if dark_shape != (2048,2048):
                
                printc("Dark Field Input File not in 2048,2048 format: {}",dark_f,color=bcolors.WARNING)
                printc("Attempting to correct ",color=bcolors.WARNING)

                
                try:
                    if dark_shape[0] > 2048:
                        dark = dark[dark_shape[0]-2048:,:]
                
453
                except Exception:
jonas's avatar
jonas committed
454
455
                    printc("ERROR, Unable to correct shape of dark field data: {}",dark_f,color=bcolors.FAIL)

Jonas Sinjan's avatar
Jonas Sinjan committed
456
457
458
            printc('--------------------------------------------------------------',bcolors.OKGREEN)
            printc(f"------------ Load darks time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
            printc('--------------------------------------------------------------',bcolors.OKGREEN)
459

jonas's avatar
jonas committed
460
461
462
        except Exception:
            printc("ERROR, Unable to open darks file: {}",dark_f,color=bcolors.FAIL)

jonas's avatar
jonas committed
463

464
465
466
        #-----------------
        # APPLY DARK CORRECTION 
        #-----------------    
Jonas Sinjan's avatar
Jonas Sinjan committed
467
        print(" ")
jonas's avatar
jonas committed
468
        print("-->>>>>>> Subtracting dark field")
Jonas Sinjan's avatar
Jonas Sinjan committed
469
        
470
471
        start_time = time.time()

Jonas Sinjan's avatar
Jonas Sinjan committed
472
        flat -= dark[..., np.newaxis, np.newaxis]
jonas's avatar
jonas committed
473

474
475
476
477
        if header_imgdirx_exists:
            if imgdirx_flipped == 'YES':
                dark_copy = np.copy(dark)
                dark_copy = dark_copy[:,::-1]
jonas's avatar
jonas committed
478

479
                data -= dark_copy[start_row:start_row + data_size[0],start_col:start_col + data_size[1], np.newaxis, np.newaxis, np.newaxis] 
jonas's avatar
jonas committed
480
481
482
            
            elif imgdirx_flipped == 'NO':
                data -= dark[start_row:start_row + data_size[0],start_col:start_col + data_size[1], np.newaxis, np.newaxis, np.newaxis] 
jonas's avatar
jonas committed
483
484
        else:
            data -= dark[start_row:start_row + data_size[0],start_col:start_col + data_size[1], np.newaxis, np.newaxis, np.newaxis] 
485

Jonas Sinjan's avatar
Jonas Sinjan committed
486
487
488
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- Dark Field correction time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
jonas's avatar
jonas committed
489

490
491
492
    else:
        print(" ")
        printc('-->>>>>>> No dark mode',color=bcolors.WARNING)
jonas's avatar
jonas committed
493

jonas's avatar
jonas committed
494

jonas's avatar
jonas committed
495
496
497
498
    #-----------------
    # OPTIONAL Unsharp Masking clean the flat field stokes V images
    #-----------------

499
    if clean_f and flat_c:
jonas's avatar
jonas committed
500
501
502
503
504
505
506
        print(" ")
        printc('-->>>>>>> Cleaning flats with Unsharp Masking (Stokes V only)',color=bcolors.OKGREEN)

        start_time = time.time()

        #cleaning the stripe in the flats for a particular flat

jonas's avatar
jonas committed
507
508
509
        # if flat_f[-62:] == 'solo_L0_phi-hrt-flat_0667134081_V202103221851C_0162201100.fits': 
        #     flat[1345, 296:, 1, 2] = flat[1344, 296:, 1, 2]
        #     flat[1346, :291, 1, 2] = flat[1345, :291, 1, 2]
jonas's avatar
jonas committed
510
511
512

        #demod the flats

jonas's avatar
jonas committed
513
        flat_demod, demodM = demod_hrt(flat, flat_pmp_temp)
jonas's avatar
jonas committed
514

515
        norm_factor = np.mean(flat_demod[512:1536,512:1536,0,0])
jonas's avatar
jonas committed
516
517
518
519
520
521
522
523
524
525
526
527
528

        flat_demod /= norm_factor

        new_demod_flats = np.copy(flat_demod)

        b_arr = np.zeros((2048,2048,3,5))

        if cpos_arr[0] == 0:
            wv_range = range(1,6)

        elif cpos_arr[0] == 5:
            wv_range = range(5)

jonas's avatar
jonas committed
529
530
531
532
533
534
535
536
        if clean_mode == "QUV":
            start_clean_pol = 1
        elif clean_mode == "UV":
            start_clean_pol = 2
        elif clean_mode == "V":
            start_clean_pol = 3

        for pol in range(start_clean_pol,4):
jonas's avatar
jonas committed
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555

            for wv in wv_range: #not the continuum

                a = np.copy(np.clip(flat_demod[:,:,pol,wv], -0.02, 0.02))
                b = a - gaussian_filter(a,sigma)
                b_arr[:,:,pol-1,wv-1] = b
                c = a - b

                new_demod_flats[:,:,pol,wv] = c

        invM = np.linalg.inv(demodM)

        flat = np.matmul(invM, new_demod_flats*norm_factor)


        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- Cleaning flat time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)

556
557
558
    else:
        print(" ")
        printc('-->>>>>>> No clean flats mode',color=bcolors.WARNING)
jonas's avatar
jonas committed
559
560


561
562
563
564
565
    #-----------------
    # NORM FLAT FIELDS
    #-----------------

    if norm_f and flat_c:
Jonas Sinjan's avatar
Jonas Sinjan committed
566
567
        
        print(" ")
568
569
570
571
572
        printc('-->>>>>>> Normalising Flats',color=bcolors.OKGREEN)

        start_time = time.time()

        try:
Jonas Sinjan's avatar
Jonas Sinjan committed
573
574
            norm_fac = np.mean(flat[512:1536,512:1536, :, :], axis = (0,1))  #mean of the central 1k x 1k
            flat /= norm_fac[np.newaxis, np.newaxis, ...]
575

Jonas Sinjan's avatar
Jonas Sinjan committed
576
577
578
            printc('--------------------------------------------------------------',bcolors.OKGREEN)
            printc(f"------------- Normalising flat time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
            printc('--------------------------------------------------------------',bcolors.OKGREEN)
579
580
581
582

        except Exception:
            printc("ERROR, Unable to normalise the flat fields: {}",flat_f,color=bcolors.FAIL)

583
584
585
586
    else:
        print(" ")
        printc('-->>>>>>> No normalising flats mode',color=bcolors.WARNING)

jonas's avatar
jonas committed
587

jonas's avatar
jonas committed
588
589
590
591
592
    #-----------------
    # APPLY FLAT CORRECTION 
    #-----------------

    if flat_c:
Jonas Sinjan's avatar
Jonas Sinjan committed
593
        print(" ")
jonas's avatar
jonas committed
594
        printc('-->>>>>>> Correcting Flatfield',color=bcolors.OKGREEN)
595
596
597

        start_time = time.time()

jonas's avatar
jonas committed
598
599
600
        if header_imgdirx_exists:
            if imgdirx_flipped == 'YES': #should be YES for any L1 data, but mistake in processing software
                flat = flat[:,::-1, :, :]
jonas's avatar
jonas committed
601
                print("Flat is flipped to match the Science Data")
jonas's avatar
jonas committed
602

jonas's avatar
jonas committed
603
        try:
604
605

            if flat_states == 6:
jonas's avatar
jonas committed
606
        
607
608
609
610
                printc("Dividing by 6 flats, one for each wavelength",color=bcolors.OKGREEN)
                    
                tmp = np.mean(flat,axis=-2) #avg over pol states for the wavelength

jonas's avatar
jonas committed
611
                data /= tmp[start_row:start_row + data_size[0],start_col:start_col + data_size[1], np.newaxis, :, np.newaxis]
612
613
614
615
616
617
618
619
620
621
622
623
624
625


            elif flat_states == 24:

                printc("Dividing by 24 flats, one for each image",color=bcolors.OKGREEN)

                data /= flat[start_row:start_row + data_size[0],start_col:start_col + data_size[1], :, :, np.newaxis] #only one new axis for the scans
                    
            elif flat_states == 4:

                printc("Dividing by 4 flats, one for each pol state",color=bcolors.OKGREEN)

                tmp = np.mean(flat,axis=-1) #avg over wavelength

jonas's avatar
jonas committed
626
                data /= tmp[start_row:start_row + data_size[0],start_col:start_col + data_size[1], :, np.newaxis, np.newaxis]
jonas's avatar
jonas committed
627
        
Jonas Sinjan's avatar
Jonas Sinjan committed
628
629
630
            printc('--------------------------------------------------------------',bcolors.OKGREEN)
            printc(f"------------- Flat Field correction time: {np.round(time.time() - start_time,3)} seconds ",bcolors.OKGREEN)
            printc('--------------------------------------------------------------',bcolors.OKGREEN)
jonas's avatar
jonas committed
631
632
633
634

        except: 
          printc("ERROR, Unable to apply flat fields",color=bcolors.FAIL)

635
636
637
638
    else:
        print(" ")
        printc('-->>>>>>> No flat field correction mode',color=bcolors.WARNING)

jonas's avatar
jonas committed
639

jonas's avatar
jonas committed
640
641
642
643
644
    #-----------------
    # PREFILTER CORRECTION  
    #-----------------

    if prefilter_f is not None:
jonas's avatar
jonas committed
645
        print(" ")
jonas's avatar
jonas committed
646
647
648
649
650
651
652
653
654
655
656
657
        printc('-->>>>>>> Prefilter Correction ',color=bcolors.OKGREEN)

        start_time = time.time()

        prefilter_voltages = [-1300.00,-1234.53,-1169.06,-1103.59,-1038.12,-972.644,-907.173,-841.702,-776.231,-710.760,-645.289,
                            -579.818,-514.347,-448.876,-383.404,-317.933,-252.462,-186.991,-121.520,-56.0490,9.42212,74.8932,
                            140.364,205.835,271.307, 336.778,402.249,467.720,533.191,598.662,664.133,729.604,795.075,860.547,
                            926.018,991.489,1056.96,1122.43,1187.90,1253.37, 1318.84,1384.32,1449.79,1515.26,1580.73,1646.20,
                            1711.67,1777.14,1842.61]

        prefilter, _ = load_fits(prefilter_f)

jonas's avatar
jonas committed
658
        #prefilter = prefilter[:,652:1419,613:1380] #crop the helioseismology data
jonas's avatar
jonas committed
659

jonas's avatar
jonas committed
660
661
662
        def get_v1_index1(x):
            index1, v1 = min(enumerate([abs(i) for i in x]), key=itemgetter(1))
            return  v1, index1
jonas's avatar
jonas committed
663

jonas's avatar
jonas committed
664
        for scan in range(data_shape[-1]):
jonas's avatar
jonas committed
665

jonas's avatar
jonas committed
666
            voltage_list = voltagesData_arr[scan]
jonas's avatar
jonas committed
667
668
            
            for wv in range(6):
jonas's avatar
jonas committed
669

jonas's avatar
jonas committed
670
                v = voltage_list[wv]
jonas's avatar
jonas committed
671

jonas's avatar
jonas committed
672
673
674
675
676
677
678
679
680
681
682
683
                vdif = [v - pf for pf in prefilter_voltages]
                
                v1, index1 = get_v1_index1(vdif)
                
                if vdif[index1] >= 0:
                    v2 = vdif[index1 + 1]
                    index2 = index1 + 1
                    
                else:
                    v2 = vdif[index1-1]
                    index2 = index1 - 1
                    
jonas's avatar
jonas committed
684
                imprefilter = (prefilter[:,:, index1]*v1 + prefilter[:,:, index2]*v2)/(v1+v2)
jonas's avatar
jonas committed
685

jonas's avatar
jonas committed
686
                data[:,:,:,wv,scan] /= imprefilter[...,np.newaxis]
jonas's avatar
jonas committed
687
688
689
690
691
692


        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- Prefilter correction time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)

693
694
695
    else:
        print(" ")
        printc('-->>>>>>> No prefilter mode',color=bcolors.WARNING)
jonas's avatar
jonas committed
696

697
698
699
700
    #-----------------
    # FIELD STOP 
    #-----------------

jonas's avatar
jonas committed
701
    if fs_c:
jonas's avatar
jonas committed
702
    
703
704
        print(" ")
        printc("-->>>>>>> Applying field stop",color=bcolors.OKGREEN)
705

706
707
708
        start_time = time.time()
        
        field_stop,_ = load_fits('./field_stop/HRT_field_stop.fits')
709

710
        field_stop = np.where(field_stop > 0,1,0)
711

jonas's avatar
jonas committed
712
713
714
        if header_imgdirx_exists:
            if imgdirx_flipped == 'YES': #should be YES for any L1 data, but mistake in processing software
                field_stop = field_stop[:,::-1] #also need to flip the flat data after dark correction
jonas's avatar
jonas committed
715

716
717
718
719
720
721
722
723
724
        data *= field_stop[start_row:start_row + data_size[0],start_col:start_col + data_size[1],np.newaxis, np.newaxis, np.newaxis]

        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- Field stop time: {np.round(time.time() - start_time,3)} seconds",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)

    else:
        print(" ")
        printc('-->>>>>>> No field stop mode',color=bcolors.WARNING)
725

jonas's avatar
jonas committed
726

jonas's avatar
jonas committed
727
728
729
730
    #-----------------
    # APPLY DEMODULATION 
    #-----------------

731
732
    if demod:

Jonas Sinjan's avatar
Jonas Sinjan committed
733
734
        print(" ")
        printc('-->>>>>>> Demodulating data         ',color=bcolors.OKGREEN)
jonas's avatar
jonas committed
735

jonas's avatar
jonas committed
736
737
        start_time = time.time()

jonas's avatar
jonas committed
738
        data,_ = demod_hrt(data, pmp_temp)
Jonas Sinjan's avatar
Jonas Sinjan committed
739
740
741
742
        
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- Demodulation time: {np.round(time.time() - start_time,3)} seconds ",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
jonas's avatar
jonas committed
743

744
745
746
    else:
        print(" ")
        printc('-->>>>>>> No demod mode',color=bcolors.WARNING)
jonas's avatar
jonas committed
747
748


jonas's avatar
jonas committed
749
750
751
752
    #-----------------
    # APPLY NORMALIZATION 
    #-----------------

753
    if norm_stokes:
Jonas Sinjan's avatar
Jonas Sinjan committed
754
755
756
757
        
        print(" ")
        printc('-->>>>>>> Normalising Stokes to Quiet Sun',color=bcolors.OKGREEN)
        
jonas's avatar
jonas committed
758
759
        start_time = time.time()

Jonas Sinjan's avatar
Jonas Sinjan committed
760
761
        for scan in range(data_shape[-1]):
            
jonas.sinjan's avatar
jonas.sinjan committed
762
            I_c = np.mean(data[512:1536,512:1536,0,cpos_arr[0],int(scan)]) #mean of central 1k x 1k of continuum stokes I
Jonas Sinjan's avatar
Jonas Sinjan committed
763
764
765
766
767
            data[:,:,:,:,scan] = data[:,:,:,:,scan]/I_c
       
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- Stokes Normalising time: {np.round(time.time() - start_time,3)} seconds ",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
jonas's avatar
jonas committed
768

769
770
771
    else:
        print(" ")
        printc('-->>>>>>> No normalising Stokes mode',color=bcolors.WARNING)
jonas's avatar
jonas committed
772
773


jonas's avatar
jonas committed
774
775
776
    #-----------------
    # CROSS-TALK CALCULATION 
    #-----------------
jonas's avatar
jonas committed
777

778
    if ItoQUV:
Jonas Sinjan's avatar
Jonas Sinjan committed
779
780
        
        print(" ")
jonas's avatar
jonas committed
781
        printc('-->>>>>>> Cross-talk correction I to Q,U,V ',color=bcolors.OKGREEN)
jonas's avatar
jonas committed
782
783
784
785
786
787
788
789
790
791
792
793

        start_time = time.time()

        before_ctalk_data = np.copy(data)

        num_of_scans = data_shape[-1]

        try:
            assert ctalk_params.shape == (2,3)
        except AssertionError:
            print("ctalk_params is not in the required (2,3) shape, please reconcile")

794
795
796
797
798
799
800
801
802
803
804

        for scan_hdr in hdr_arr:
            if 'CAL_CRT0' in scan_hdr: #check to make sure the keywords exist
                scan_hdr['CAL_CRT0'] = ctalk_params[0,0] #I-Q slope
                scan_hdr['CAL_CRT2'] = ctalk_params[0,1] #I-U slope
                scan_hdr['CAL_CRT4'] = ctalk_params[0,2] #I-V slope

                scan_hdr['CAL_CRT1'] = ctalk_params[1,0] #I-Q offset
                scan_hdr['CAL_CRT3'] = ctalk_params[1,1] #I-U offset
                scan_hdr['CAL_CRT5'] = ctalk_params[1,2] #I-V offset

jonas's avatar
jonas committed
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
        ctalk_params = np.repeat(ctalk_params[:,:,np.newaxis], num_of_scans, axis = 2)

        cont_stokes = np.mean(data[512:1536,512:1536,0,cpos_arr[0],:], axis = (0,1))

        for i in range(6):
            stokes_i_wv_avg = np.mean(data[512:1536,512:1536,0,i,:], axis = (0,1))

            if norm_stokes:
                #if normed, applies normalised offset to normed stokes

                tmp_param = ctalk_params*np.divide(stokes_i_wv_avg,cont_stokes)

                q_slope = tmp_param[0,0,:]
                u_slope = tmp_param[0,1,:]
                v_slope = tmp_param[0,2,:]

                q_int = tmp_param[1,0,:]
                u_int = tmp_param[1,1,:]
                v_int = tmp_param[1,2,:]

                data[:,:,1,i,:] = before_ctalk_data[:,:,1,i,:] - before_ctalk_data[:,:,0,i,:]*q_slope - q_int

                data[:,:,2,i,:] = before_ctalk_data[:,:,2,i,:] - before_ctalk_data[:,:,0,i,:]*u_slope - u_int

                data[:,:,3,i,:] = before_ctalk_data[:,:,3,i,:] - before_ctalk_data[:,:,0,i,:]*v_slope - v_int

            else:
                #if not normed, applies raw offset cross talk correction to raw stokes counts

                tmp_param = ctalk_params[0,:,:]*np.divide(stokes_i_wv_avg,cont_stokes)

                q_slope = tmp_param[0,:]
                u_slope = tmp_param[1,:]
                v_slope = tmp_param[2,:]

                q_int = ctalk_params[1,0,:]
                u_int = ctalk_params[1,1,:]
                v_int = ctalk_params[1,2,:]

                data[:,:,1,i,:] = before_ctalk_data[:,:,1,i,:] - before_ctalk_data[:,:,0,i,:]*q_slope - q_int*stokes_i_wv_avg 

                data[:,:,2,i,:] = before_ctalk_data[:,:,2,i,:] - before_ctalk_data[:,:,0,i,:]*u_slope - u_int*stokes_i_wv_avg 

                data[:,:,3,i,:] = before_ctalk_data[:,:,3,i,:] - before_ctalk_data[:,:,0,i,:]*v_slope - v_int*stokes_i_wv_avg 


jonas's avatar
jonas committed
851
852
853
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
        printc(f"------------- I -> Q,U,V cross talk correction time: {np.round(time.time() - start_time,3)} seconds ",bcolors.OKGREEN)
        printc('--------------------------------------------------------------',bcolors.OKGREEN)
jonas's avatar
jonas committed
854
855
856
        
        data *= field_stop[start_row:start_row + data_size[0],start_col:start_col + data_size[1], np.newaxis, np.newaxis, np.newaxis]

857
858
859
860
    else:
        print(" ")
        printc('-->>>>>>> No ItoQUV mode',color=bcolors.WARNING)

jonas's avatar
jonas committed
861
862
863
864
865
866
867
868

    #-----------------
    #CHECK FOR INFs
    #-----------------

    data[np.isinf(data)] = 0
    data[np.isnan(data)] = 0

jonas.sinjan's avatar
jonas.sinjan committed
869
    
870
    if out_demod_file:
871
872
873

        #check if the output directory exists, if not, create it
        if not os.path.exists(out_dir): 
jonas's avatar
jonas committed
874
            print(f"{out_dir} does not exist -->>>>>>> Creating it")
875
            os.makedirs(out_dir)
876
        
877
878
        print(" ")
        printc('Saving demodulated data to one _reduced.fits file per scan')
879

880
        if out_demod_filename is not None:
jonas's avatar
jonas committed
881
882
883
884
885

            if isinstance(out_demod_filename,str):
                out_demod_filename = [out_demod_filename]

            if int(len(out_demod_filename)) == int(data_shape[-1]):
jonas's avatar
jonas committed
886
                scan_name_list = out_demod_filename
887
888
889
890
891
892
893
894
895
896
                scan_name_defined = True
            else:
                print("Input demod filenames do not match the number of input arrays, reverting to default naming")
                scan_name_defined = False
        else:
            scan_name_defined = False

        if not scan_name_defined: #check if already defined by input, otherwise generate
            scan_name_list = check_filenames(data_f)
        
jonas's avatar
jonas committed
897

898
        for count, scan in enumerate(data_f):
899

900
            with fits.open(scan) as hdu_list:
jonas's avatar
jonas committed
901
                print(f"Writing out demod file as: {scan_name_list[count]}_reduced.fits")
902
                hdu_list[0].data = data[:,:,:,:,count]
903
                hdu_list[0].header = hdr_arr[count] #update the calibration keywords
904
905
906
907
908
909
910
911
912
913
914
                hdu_list.writeto(out_dir + scan_name_list[count] + '_reduced.fits', overwrite=True)

        # if isinstance(data_f, str):
        #     print(" ")
        #     printc('Saving demodulated data to a _reduced.fits file')

        #     if out_demod_filename is not None:
        #         scan_name = str(out_demod_filename)
            
        #     if 'scan_name' not in locals(): #check if already defined by input, otherwise generate
        #         scan_name = str(data_f.split('.fits')[0][-10:])
915

916
        #     with fits.open(data_f) as hdu_list:
917

918
919
        #         hdu_list[0].data = data
        #         hdu_list.writeto(out_dir + scan_name + '_reduced.fits', overwrite=True)
jonas's avatar
jonas committed
920

921
922
    else:
        print(" ")
923
924
        #check if already defined by input, otherwise generate
        scan_name_list = check_filenames(data_f)
925
        printc('-->>>>>>> No output demod file mode',color=bcolors.WARNING)
jonas's avatar
reorg    
jonas committed
926

jonas's avatar
jonas committed
927
928
929
930
931
    #-----------------
    # INVERSION OF DATA WITH CMILOS
    #-----------------

    if rte == 'RTE' or rte == 'CE' or rte == 'CE+RTE':
jonas's avatar
jonas committed
932

933
        if p_milos:
jonas's avatar
jonas committed
934

jonas's avatar
jonas committed
935
936
937
938
939
            try:
                pmilos(data_f, wve_axis_arr, data_shape, cpos_arr, data, rte, field_stop, start_row, start_col, out_rte_filename, out_dir)
                    
            except ValueError:
                print("Running CMILOS instead!")
jonas's avatar
jonas committed
940
                cmilos(data_f, hdr_arr, wve_axis_arr, data_shape, cpos_arr, data, rte, field_stop, start_row, start_col, out_rte_filename, out_dir)
jonas's avatar
jonas committed
941

jonas's avatar
jonas committed
942
        else:
jonas's avatar
jonas committed
943
            cmilos(data_f, hdr_arr, wve_axis_arr, data_shape, cpos_arr, data, rte, field_stop, start_row, start_col, out_rte_filename, out_dir)
jonas's avatar
jonas committed
944

945
946
947
948
    else:
        print(" ")
        printc('-->>>>>>> No RTE Inversion mode',color=bcolors.WARNING)

jonas's avatar
jonas committed
949
950
951
952
953
954
955
956
957
    
    #-----------------
    # CONFIG FILE
    #-----------------

    if config_file:
        print(" ")
        printc('-->>>>>>> Writing out config file ',color=bcolors.OKGREEN)

jonas's avatar
jonas committed
958
        json.dump(saved_args, open(f"{out_dir + scan_name_list[count]}" + '_hrt_pipeline_config.txt', "w"))
jonas's avatar
jonas committed
959

Jonas Sinjan's avatar
Jonas Sinjan committed
960
961
962
963
    print(" ")
    printc('--------------------------------------------------------------',color=bcolors.OKGREEN)
    printc(f'------------ Reduction Complete: {np.round(time.time() - overall_time,3)} seconds',color=bcolors.OKGREEN)
    printc('--------------------------------------------------------------',color=bcolors.OKGREEN)
jonas's avatar
jonas committed
964
965


jonas's avatar
jonas committed
966
    return data