#!/usr/bin/env python3
import os, sys, math
sys.path.append('/auto')
from   pet_misc  import *

source = "J1215+3448"
z = 0.851

if ( len(sys.argv)-1 < 3 ):
     print ( "Usage: gen_output_model exp_date band niter" )
     exit  ( 1 )
else:
     exp_dir  = sys.argv[1]
     band     = sys.argv[2]
     niter    = int(sys.argv[3])

uvs_fits_file = exp_dir + "/" + source + "_" + band + "_" + "uvs.fits"
fit_mod_file  = exp_dir + "/" + source + "_" + band + "_" + "fit.mod"
fit1_mod_file  = exp_dir + "/" + source + "_" + band + "_" + "fit1.mod"

buf_fit = read_file ( fit_mod_file )
for line in buf_fit:
    if ( line[0:1] != "!" ):
         buf_fit1 = [ "%s 0.0 0.0 2.0v 1.0v 0.0v 1" % line.split()[0] ]
         break

(ret,err) = write_file ( buf_fit1, fit1_mod_file )

obs_date = exp_dir.split("_")[0]
exp_name = exp_dir.split("_")[1]

if ( not os.path.isfile ( uvs_fits_file ) ):
     print ( "Cannot find calibrated visibility file %s" % uvs_fits_file )
     exit  ( 1 )

if ( not os.path.isfile ( fit_mod_file ) ):
     print ( "Cannot find model file %s" % fit_mod_file )
     exit  ( 1 )

cnt_file_buf = []
cnt_file_buf.append ( "obs  %s" % uvs_fits_file )
cnt_file_buf.append ( "rmod %s" % fit_mod_file  )
cnt_file_buf.append ( "sele I"                  )
cnt_file_buf.append ( "modelfit %2d" % niter    )
cnt_file_buf.append ( "mapsize 1024,0.4, 1024,0.4;clean 1,0.05,0.0"        )
cnt_file_buf.append ( "quit"                    )

cnt_file = "/tmp/gen_mod__%08d.cnt" % os.getpid()

(ret,err) = write_file ( cnt_file_buf, cnt_file )
check_err_exe ( ret, err, "write_file" )

com = "/opt64/bin/difmap < " + cnt_file

(ret,out) = exe ( com )
os.unlink  ( cnt_file )

if ( ret != 0 ):
     for line in out:
         print ( line )
     print ( "ERROR in executing command %s" % com )
     exit  ( 1 )

cnt_file_buf = []
cnt_file_buf.append ( "obs  %s" % uvs_fits_file )
cnt_file_buf.append ( "rmod %s" % fit1_mod_file )
cnt_file_buf.append ( "sele I"                  )
cnt_file_buf.append ( "modelfit %2d" % niter    )
cnt_file_buf.append ( "mapsize 1024,0.4, 1024,0.4;clean 1,0.05,0.0"        )
cnt_file_buf.append ( "quit"                    )

(ret,err) = write_file ( cnt_file_buf, cnt_file )
check_err_exe ( ret, err, "write_file" )

(ret,out1) = exe ( com )
# os.unlink  ( cnt_file )

if ( ret != 0 ):
     for line in out1:
         print ( line )
     print ( "ERROR in executing command %s" % com )
     exit  ( 1 )

res = []
for i in range(0,len(out)):
    if ( " gauss " in out[i] ):
         flux_B     = float ( out[i].split()[0] )
         flux_B_err = float ( out[i].split()[1] )
         flux_B_err = max ( flux_B_err, 0.05*flux_B )
         dx         = float ( out[i].split()[2] ) *1000.0
         dx_e       = float ( out[i].split()[3] ) *1000.0
         dy         = float ( out[i].split()[4] ) *1000.0
         dy_e       = float ( out[i].split()[5] ) *1000.0
         dr         = math.sqrt ( dx**2 + dy**2 )
         dr_e       = (abs(dx)*dx_e + abs(dy)*dy_e)/dr
         if ( exp_name == "bb023g" or exp_name == "rdv36" ):
              dx_e = 0.001*dx_e
              dy_e = 0.001*dy_e
              dr_e = 0.001*dr_e
              flux_B_err = 0.001*flux_B_err
         dr_e = max ( dr_e, 0.01 )

    if ( "Estimated beam:" in out[i] ):
          beam_min = float ( out[i].split()[2].replace("bmin=","" ) )
          beam_maj = float ( out[i].split()[4].replace("bmaj=","" ) )
          beam_pa  = float ( out[i].split()[6].replace("bpa=","" ) )

for i in range(0,len(out1)):
    line = out1[i]
    if ( "Quitting program" in out1[i] ):
         line = out1[i].replace("Quitting program"," " ) + out1[i+1]
    if ( " gauss " in out1[i] ):
          flux_A_peak     = float ( line.split()[0] )
          flux_A_peak_err = float ( line.split()[1] )
          if ( exp_name == "bb023g" or exp_name == "rdv36" ):
               flux_A_peak_err = 0.001*flux_A_peak_err  # because of AIPS
          flux_A_peak_err = max ( flux_A_peak_err, 0.05*flux_A_peak )
          A_peak_min  = 1000.0 * float ( line.split()[9].replace("bmin=","" )  )
          A_peak_maj  = 1000.0 * float ( line.split()[11].replace("bmaj=","" ) )
          A_peak_bpa  =          float ( line.split()[13].replace("bpa=","" )  )
          break

mod_buf = read_file ( fit_mod_file )
for i in range(0,len(mod_buf)):
    if ( "Total flux:" in mod_buf[i] ):
         tot_flux = float ( mod_buf[i].split()[3] )

# print ( "flux_A_peak_err ", flux_A_peak_err ) # %%%%%%%%%%%%%%

if ( band == "L" ): 
     freq = 1.6
elif ( band == "S" ): 
     freq = 2.3
elif ( band == "C" ): 
     freq = 4.3
     if ( exp_name == "bh019" or exp_name == "bv019b" ):
          freq = 5.0 
elif ( band == "X" ): 
     freq = 8.4
     if ( exp_name == "bp175a1" or exp_name == "bm425a" or exp_name == "sb072d1" ):
          freq = 7.6
elif ( band == "U" ): 
     freq = 15.2
elif ( band == "K" ): 
     freq = 22.2
elif ( band == "Q" ): 
     freq = 43.0


snr = flux_A_peak / flux_A_peak_err
A_lim_min = beam_min * math.sqrt ( 4.0*math.log(2.0)/3.14159 * math.log(snr/(snr-1)) )
A_lim_maj = beam_maj * math.sqrt ( 4.0*math.log(2.0)/3.14159 * math.log(snr/(snr-1)) )
peak_a_area = max ( A_peak_maj, A_lim_maj ) * max ( A_peak_min, A_lim_min )
if ( A_peak_maj < A_lim_maj or A_peak_min < A_lim_min ):
     sign_tb_a = ">"
else:
     sign_tb_a = " "

Tb_A = 1.22E12 * flux_A_peak/(freq**2 * peak_a_area)*(1.0 + z)
Tb_B = 1.22E12 * flux_B/(freq**2 * beam_min*beam_maj/4)*(1.0 + z)

print ( "Obs_date: %s  Exp_name: %-7s  Band: %s Total_flux: %6.3f  Tot_flux_A: %6.3f Flux_A_peak: %6.3f Flux_B: %6.3f -+ %6.3f Jy  Dx: %8.2f -+ %5.2f mas  Dy: %8.2f -+ %5.2f mas  Dist: %6.2f -+ %5.2f mas  Beam: %4.1f x %4.1f bpa: %6.1f A_unr_min: %4.1f A_unr_maj: %4.1f A_bpa: %6.1f &  Tb_A: %s%.1E  Tb_B: >%.1E" % \
        ( obs_date, exp_name, band, tot_flux, tot_flux-flux_B, flux_A_peak, flux_B, flux_B_err, \
          dx, max(dx_e,0.01), dy, max(dy_e,0.01), dr, dr_e, max(beam_min,0.1), max(beam_maj,0.1), beam_pa, max(A_peak_min,0.1), max(A_peak_maj,0.1), A_peak_bpa, sign_tb_a, Tb_A, Tb_B ) )

print ( "TEX %s & %-7s & %4.1f & %5.3f & %6.3f & %6.3f & $\\pm$ %5.3f & %5.3f & $\\pm$ %5.3f  %6.2f & $\\pm$ %4.2f & %6.2f & $\\pm$ %4.2f & %4.1f $\\times$ %4.1f & %6.1f & %4.1f $\\times$  %4.1f & %6.1f & $%s$%.1E  & $>$%.1E \\\\" % \
        ( obs_date, exp_name, freq, tot_flux, tot_flux-flux_B, flux_A_peak, flux_A_peak_err, flux_B, flux_B_err, \
          dx, max(dx_e,0.01), dy, max(dy_e,0.01), max(beam_min,0.1), max(beam_maj,0.1), beam_pa, max(A_peak_min,0.1), max(A_peak_maj,0.1), A_peak_bpa, sign_tb_a, Tb_A, Tb_B ) )