#!/usr/bin/perl -w
use strict;
# variable
	my($ra1,$dec1,$ra2,$dec2);
	my(@list,@t,$src);
	my($r1,$r2,$r3,$a1,$a2,$a3,$dia,$nu,$d,$f);
	my($r11,$r22,$r33);
	my($nvss1,$nvss2,$radius1,$radius2,$flux1,$flux2,$f1,$f2);
	my($range1,$range2);
	my($rms,$sigma,$band,$dynamic_range1,$dynamic_range2,$obs_type);

# inputs
	# printf(" * $ARGV[0] $ARGV[1] $ARGV[2] $ARGV[3] $ARGV[4] $ARGV[5] $ARGV[6]\n"); exit 0;

# inputs from php
	$ra1=$ARGV[0]; # radian
	$dec1=$ARGV[1]; # radian

	$dia=$ARGV[2]; # meter
	$nu=$ARGV[3]; # GHz

	$rms=$ARGV[4]; # microJy
	$sigma=$ARGV[5]; # microJy
	$band=$ARGV[6];  # band (Band2, Band3, Band4, Band5)
	$obs_type=$ARGV[7];

# calculate r1, r2, r3
	# r1
	$a1=37102.06877226444;
	$r1=($a1/($nu*$dia));  # arc sec
	$r11=sprintf("%d",$r1/60); # arc min
	$r1=($r1/3600)/57.29577951308232;  # radian 

	# r2
	$a2=35250.0;
	$r2=($a2/($nu*$dia));  # arc sec
	$r22=sprintf("%d",$r2/60); # arc min
	$r2=($r2/3600)/57.29577951308232;  # radian 

	# r3
	$a3=44522.0;
	$r3=($a3/($nu*$dia));  # arc sec
	$r33=sprintf("%d",$r3/60); # arc min
	$r3=($r3/3600)/57.29577951308232;  # radian 

# NVSS Dec search range (+/- 3deg)
	$range1=$dec1-0.0523599; # radian
	$range2=$dec1+0.0523599; # radian

# read nvss list
	chomp(@list=`cat nvss_radian.dat | awk ' { if( (\$3 > $range1) && (\$3 < $range2)) { print \$0; }}' `);

# search in the list 
	$nvss1=0;
	$nvss2=0;
	$radius1=-1;
	$radius2=-1;
	$flux1=0;
	$flux2=0;
	$f1=0;
	$f2=0;

foreach $src (@list) { 
	@t=split(/ +/,$src);

	$f=$t[0];
	$ra2=$t[1];
	$dec2=$t[2];

	# calculate the angular distance between target and caltalog source 
	$d=sin($dec1)*sin($dec2)+cos($dec1)*cos($dec2)*cos($ra1-$ra2);
	$d=acos($d); 

	# search between 0 - r1
	if((!$nvss1) && ($d < $r1) && ($f > 50)) { 
	$nvss1=$src;
	$radius1=sprintf("%d",$d*57.29577951308232*60); # arcmin
	$flux1=sprintf("%2.1f",$f*(($nu/1.4)**(-0.7)));
	$f1=$f;
	}

	# search between r2 - r3
	if((!$nvss2) && ($d >= $r2) && ($d <= $r3) && ($f > 100)) { 
	$nvss2=$src;
	$radius2=sprintf("%d",$d*57.29577951308232*60); # arcmin
	$flux2=sprintf("%2.1f",$f*(($nu/1.4)**(-0.7)));
	$f2=$f;
	}

	# end search
	if($nvss1 && $nvss2) { last; } 
}

# calculate the image dynamic range for source1
	if($nvss1){ 
	if($band eq "Band3" || $band eq "Band4" || $band eq "Band5" ){
	$dynamic_range1=$flux1*1000/$rms;
	}
	if($band eq "Band2"){
	if($rms >= $sigma) { $dynamic_range1=$flux1*1000/$rms; } 
	if($rms <  $sigma) { $dynamic_range1=$flux1*1000/$sigma; } 
	}
	$dynamic_range1=sprintf("%d",$dynamic_range1);
	}

# calculate the image dynamic range for source2
	if($nvss2){ 
	if($band eq "Band3" || $band eq "Band4" || $band eq "Band5" ){
	$dynamic_range2=$flux1*1000/$rms;
	}
	if($band eq "Band2"){
	if($rms >= $sigma) { $dynamic_range2=$flux2*1000/$rms; } 
	if($rms <  $sigma) { $dynamic_range2=$flux2*1000/$sigma; } 
	}
	$dynamic_range2=sprintf("%d",$dynamic_range2);
	}

# message1
	if(!$nvss1 && !$nvss2) { 
	printf("\n");
	printf("There is no source brighter than 50 mJy in the 1.4 GHz NVSS within the FWHM, of the GMRT primary beam at the observing frequency."); 
	printf("\n");
	printf("There is also no source brighter than 100 mJy between the FWHM and the first null of the primary beam at the observing frequency.");
	}

# message2
	if($nvss1 && !$nvss2) { 
	@t=split(/ +/,$nvss1);
	printf("\n");
	printf("The brightest source within the FWHM of the GMRT primary beam at the observing frequency has an extrapolated flux density of $flux1 mJy.");
	printf(" It is located at RA = ($t[3] $t[4] $t[5]) and Dec = ($t[6] $t[7] $t[8]), at distance of $radius1' from the target."); 
	if($obs_type =~ "Continuum") { 
	printf(" The dynamic range of the image would be ~ $dynamic_range1."); 
	}
	printf("\n");
	printf(" There is no source brighter than 100 mJy in the 1.4 GHz NVSS"); 
	printf(" between the FWHM and the first null of the GMRT primary beam at the observing frequency.");
	}

# message3
	if(!$nvss1 && $nvss2) { 
	@t=split(/ +/,$nvss2);
	printf("\n");
	printf("There is no source brighter than 50 mJy in the 1.4 GHz NVSS within the FWHM of the GMRT primary beam at the observing frequency.");
	printf("\n");
	printf(" The brightest source between the FWHM and the first null of the GMRT primary beam at the observing frequency has an extrapolated flux density of $flux2 mJy.");
	printf(" It is located at RA = ($t[3] $t[4] $t[5]) and Dec = ($t[6] $t[7] $t[8]), at a distance of $radius2' from the target.");
	printf(" The image quality may be limited by deconvolution errors.");
	}

# message4
	if($nvss1 && $nvss2) { 
	@t=split(/ +/,$nvss1);
	printf("\n");
	printf("The brightest source within the FWHM of the GMRT primary beam at the observing frequency has an extrapolated flux density of $flux1 mJy.");
	printf(" It is located at RA = ($t[3] $t[4] $t[5]) and Dec = ($t[6] $t[7] $t[8]), at a distance of $radius1' from the target.");
	if($obs_type =~ "Continuum") { 
	printf(" The dynamic range of the image would be ~ $dynamic_range1.");
	}
	printf("\n");
	@t=split(/ +/,$nvss2);
	printf("The brightest source between the FWHM and the first null of the GMRT primary beam at the observing frequency has an extrapolated flux density of $flux2 mJy.");
	printf(" It is located at RA = ($t[3] $t[4] $t[5]) and Dec = ($t[6] $t[7] $t[8]), at a distance of $radius2' from the target."); 
	printf(" The image quality may be limited by deconvolution errors.");
	}

# warning
	if($band eq "Band2"){ 
	if($dynamic_range1 > 10000) { 
	printf("\n");
	printf("WARNING:");
	printf(" The target image dynamic range is $dynamic_range1. The image quality may be limited by dynamic range issues."); 
	} }

	if($band eq "Band3"){ 
	if($dynamic_range1 > 30000) { 
	printf("\n");
	printf("WARNING:");
	printf(" The target image dynamic range is $dynamic_range1. The image quality may be limited by dynamic range issues."); 
	} }

	if($band eq "Band4" || $band eq "Band5"){ 
	if($dynamic_range1 > 50000) { 
	printf("\n");
	printf("WARNING:");
	printf(" The target image dynamic range is $dynamic_range1. The image quality may be limited by dynamic range issues."); 
	} }


# end

sub ratorad{
	use strict;
	my($ra)=@_;
	my(@t)=split(/([hms])/,$ra);
	my($deg)=$t[0]*3600.0+$t[2]*60.0+$t[4];
	$deg=$deg/240.0;
	my($rad)=$deg/57.29577951308232;
	return($rad);
	}

sub dectorad{
	use strict;
        my($dec)=@_;
        my($deg,$rad,$sign,@t);
        @t=split(//,$dec);
        $sign="+";
        if($t[0] eq "-") { $sign="-";}
        @t=split(/([d'"])/,$dec);
        if($sign eq "+"){$deg=$t[0]+$t[2]/60.0+$t[4]/3600;}
        if($sign eq "-"){$deg=$t[0]-$t[2]/60.0-$t[4]/3600;}
        $rad=$deg/57.29577951308232;
        return($rad);
        }


sub radtodeg{ 
	use strict;
	my($rad) = @_;
	my($deg)=$rad*57.29577951308232;
	return $deg;
	}

sub acos { atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ) }
sub asin { atan2($_[0], sqrt(1 - $_[0] * $_[0])) }
