1 #!/usr/bin/perl
        2
        3 use PDL;
        4 use PDL::Transform;
        5 use PDL::NiceSlice;
        6 use PDL::Image2D;
        7 use PDL::Transform;
        8 use PDL::IO::Dumper;
        9 use PDL::Minuit;
       10 use PDL::GSLSF::LEGENDRE;
       11 use Astro::Coord;
       12 use Astro::Time;
       13 use TrIm;
       14 use Vtools;
       15 use strict;
       16 use warnings;
       17
       18 my $pi = 2*acos(0);
       19 #globals;
       20
       21 #my $data_pass;
       22 #my $im_pass;
       23
       24 my $h_pass;
       25
       26 #my $dumpit = 0;
       27 #my $debug = 0;
       28 #
       29 #my $offset;
       30 #my $scl;
       31 #
       32 #my $nord_x_pass;
       33 #my $nord_y_pass;
       34 #
       35 my $min_x;
       36 my $max_x;
       37 my $min_y;
       38 my $max_y;
       39 #
       40 {
       41
       42 # INPUT IMAGE
       43     my $im = rfits('/canopus/goto/110907/NGC6302_G_5s_dark.fit');
       44     my $h = $im -> gethdr;
       45 # GOTO field of view
       46     my $field_0 = 9./60.; #field in deg
       47 # GOTO pixel scale
       48     my $pixscale = $field_0 / $im->getdim(0);
       49
       50
       51 # GIVE NAMES, REF RA DEC and guess pixel location
       52     my @data = (
       53                 {Name => 'e1', RA  => '17:13:53.950', DEC => '-37:08:37.87', i_init => 1112, j_init => 44},
       54                 {Name => 'e2', RA  => '17:13:55.939', DEC => '-37:08:05.55', i_init => 1204, j_init => 141},
       55                 {Name => 'e3', RA  => '17:13:52.242', DEC => '-37:07:34.97', i_init => 1071, j_init => 233},
       56                 {Name => 'e4', RA  => '17:13:53.947', DEC => '-37:07:17.97', i_init => 1137, j_init => 279},
       57                 {Name => 'e5', RA  => '17:13:57.928', DEC => '-37:07:23.04', i_init => 1270, j_init => 258},
       58                 {Name => 'e6', RA  => '17:14:00.485', DEC => '-37:07:09.42', i_init =>1362 , j_init => 289},
       59                 {Name => 'e7', RA  => '17:14:01.483', DEC => '-37:07:48.51', i_init => 1388, j_init => 166},
       60                 {Name => 'e8', RA  => '17:13:52.236', DEC => '-37:04:09.27', i_init => 1122, j_init => 831},
       61                 {Name => 'e9', RA  => '17:13:31.925', DEC => '-37:04:00.75', i_init => 416, j_init => 908},
       62
       63                 {Name => 'e10', RA  => '17:13:33.059', DEC => '-37:05:10.46', i_init => 432, j_init => 703},
       64                 {Name => 'e11', RA  => '17:13:26.098', DEC => '-37:05:05.30', i_init => 196, j_init => 734},
       65                 {Name => 'e13', RA  => '17:13:24.252', DEC => '-37:04:43.19', i_init => 135, j_init => 800},
       66                 {Name => 'e14', RA  => '17:13:25.382', DEC => '-37:06:16.70', i_init => 155, j_init => 519},
       67                 );
       68
       69
       70
       71 # Define a center for the resampled image
       72
       73     my $RA_0 = '17:13:44.282'; my $DEC_0 = '-37:06:15.09'; #NGC 4755
       74
       75 # Define  dimensions for the resamepld image
       76     my $Nx = 1800;
       77     my $Ny = 1300;
       78
       79     my $nord_x = 2;
       80     my $nord_y = 1;
       81 # get resampled image:
       82     my $out = &astrogoto($im,{
       83         Nx => $Nx,
       84         Ny => $Ny,
       85         PixScale => $pixscale,
       86         RA_0 => $RA_0,
       87         DEC_0=>  $DEC_0,
       88         Nord_x => $nord_x,
       89         Nord_y => $nord_y,
       90         Data => \@data
       91     });
       92
       93     wfits $out,'testG.fits'; # choose an ouput filename
       94
       95 }
       96
       97
       98
       99
      100 sub get_centroid {
      101     my ($im, $x0, $y0, $opts) = @_;
      102
      103     my %args = ( Do_Gaussfit => 0, Delta => 6, %$opts);
      104
      105 #    if (ref($opts) =~ /HASH/) {
      106 #       foreach (keys(%$opts)) {
      107 #           $args{$_} = $$opts{$_};
      108 #       }
      109 #    }
      110
      111     my $delta=$args{Delta};
      112
      113     my $subim = $im($x0-$delta:$x0+$delta,$y0-$delta:$y0+$delta);
      114     my ($max,$xcent,$ycent) = max2d_ind($subim);
      115     $xcent += $x0-$delta;
      116     $ycent += $y0-$delta;
      117
      118     if ($args{Do_Gaussfit}) {
      119         my $civs = pdl [ 0.4481183, 0.62394083, 0.014348385];
      120         my $params = { Cinv => $civs ,
      121                        Flux=> $max->sclr,
      122                        Centroid=>pdl [ $xcent->sclr, $ycent->sclr]} ;
      123
      124         # ....
      125     }
      126     return ($xcent,$ycent);
      127 }
      128
      129
      130
      131
      132
      133
      134
      135 sub astrogoto {
      136     my ($im,$pass) = @_;
      137     my $h = $im-> gethdr;
      138     my %args = (
      139         Nord_x=> 1,
      140         Nord_y => 1,
      141         Nx => 0,
      142         Ny => 0 ,
      143         Data => 0,
      144         RA_0 => ,
      145         DEC_0 => ,
      146         PixScale => 0
      147     );
      148
      149     if (defined($pass)) {
      150         %args = (%args,  %$pass);
      151     }
      152     my $pixscale = $args{PixScale};
      153
      154
      155 # GIVE NAMES, REF RA DEC and guess pixel location
      156     my $datal = $args{Data};
      157     my @data = @$datal;
      158     my $RA_0 = $args{RA_0}; my $DEC_0 = $args{DEC_0};
      159
      160 # Define  dimensions for the resamepld image
      161     my $Nx = $args{Nx};
      162     my $Ny = $args{Ny};
      163
      164
      165
      166
      167     $$h{'CRVAL1'} = str2deg($RA_0,'H');
      168     $$h{'CRVAL2'} = str2deg($DEC_0,'D');
      169     $$h{'CRPIX1'} = $Nx/2;
      170     $$h{'CRPIX2'} = $Ny/2;
      171
      172     $$h{'CDELT1'} = -1 * $pixscale;
      173     $$h{'CDELT2'} = $pixscale;
      174     $$h{'CTYPE1'} = 'RA---TAN';
      175     $$h{'CTYPE2'} = 'DEC--TAN';
      176
      177 #    my $out = zeroes($im->dims);
      178     my $out = zeroes($Nx,$Ny);
      179     $$h{NAXIS1} = $Nx;
      180     $$h{NAXIS2} = $Ny;
      181     my ($Nx_in,$Ny_in) = $im->dims;
      182     $out(:$Nx_in-1,:$Ny_in-1) .= $im(:,:);
      183     $out ->sethdr($h);
      184     $h_pass = $h;
      185
      186 #    my $t_out = t_fits($out);
      187
      188     my @i_in;
      189     my @i_out;
      190
      191     my @j_in;
      192     my @j_out;
      193
      194     foreach (@data) {
      195         my ($i,$j)
      196             = &get_centroid($im,$_->{i_init},$_->{j_init}, {Delta => 20});
      197         $_->{i} = $i;
      198         $_->{j} = $j;
      199
      200 #    my $in = pdl (str2deg($_->{RA},'H'), str2deg($_->{DEC},'D'));
      201 #    my ($i_out, $j_out) = list ($in -> invert($t_out));
      202
      203         my $RA_OFF
      204             = str2deg($RA_0,'H')
      205             + (
      206                 str2deg($_->{RA},'H')
      207               -
      208                 str2deg($RA_0,'H')
      209               )
      210               * cos($pi*str2deg($DEC_0,'D') / 180.)
      211             ;
      212
      213         my ($i_out, $j_out)
      214             = &Ftools::adij_linear($out,$RA_OFF, str2deg($_->{DEC},'D'));
      215
      216         $_->{i_out} = $i_out;
      217         $_->{j_out} = $j_out;
      218
      219
      220 #    print "$i $j ---> $i_out $j_out \n";
      221         push @i_in,$i;
      222         push @i_out,$i_out;
      223
      224         push @j_in,$j;
      225         push @j_out,$j_out;
      226
      227
      228     }
      229
      230     my @dum = @i_out;
      231
      232     push @dum, (0, $$h{NAXIS1});
      233     $min_x = minimum(pdl(@dum));
      234     $max_x = maximum(pdl(@dum));
      235
      236     @dum = @j_out;
      237     push @dum, (0, $$h{NAXIS2});
      238     $min_y = minimum(pdl(@dum));
      239     $max_y = maximum(pdl(@dum));
      240
      241
      242
      243     my $i_in = pdl (@i_in);
      244     my $i_out = pdl (@i_out);
      245     wcols $i_in, $i_out;
      246
      247     print "i_out $i_out \n";
      248
      249 #    &Vtools::spec($i_in, $i_out);
      250 #    &Vtools::spec(pdl(@j_in), pdl(@j_out));
      251
      252
      253 # need to specify total range of pixel coordinates to pass argumens
      254 # as cos(theta) in Pl(cos(theta)) - may not be present in the
      255 # fits header of a single position vector
      256     my $opts
      257         = {
      258             Min_pixcoord_x => $min_x,
      259             Max_pixcoord_x => $max_x,
      260             Min_pixcoord_y => $min_y,
      261             Max_pixcoord_y => $max_y
      262         };
      263     $opts->{Nord_x} = $pass->{Nord_x};
      264     $opts->{Nord_y} = $pass->{Nord_y};
      265
      266     my $bestfit = &TrIm::opt_2DTransform(\@data, $opts);
      267
      268     foreach (@data) {
      269         my $in = pdl ( $_->{i}, $_->{j});
      270         my $out = pdl ( $_->{i_out}, $_->{j_out});
      271         my $model
      272             = &TrIm::Do_Transform(
      273                 $out,
      274                 $bestfit,
      275                 {
      276                     Min_pixcoord_x => $min_x,
      277                     Max_pixcoord_x => $max_x,
      278                     Min_pixcoord_y => $min_y,
      279                     Max_pixcoord_y => $max_y,
      280                     Dir => -1
      281                 }
      282             );
      283         print "$_->{Name} in $in  -> out $out -> model $model CRPIX1 $$h{CRPIX1} CRPIX2 $$h{CRPIX2} \n";
      284
      285 #       my $inversemodel = &Do_Transform($model,$bestfit,{Dir => -1});
      286 #       print "-> inversemodel $inversemodel    \n";
      287     }
      288
      289 #    my $testin = $data[-1];
      290 #    my $t_in_0 = pdl ( $testin->{i}, $testin->{j});
      291 #
      292 #    my $testout = &Do_Transform($t_in_0,$bestfit,{Dir => 1, Debug => 1});
      293 #
      294 #    print "testin $testin->{Name} $t_in_0  ---> testout $testout   \n";
      295 #
      296 #    my $t_in = pdl ( $testin->{i_out}, $testin->{j_out});
      297 #
      298 #    $testout = &Do_Transform($t_in,$bestfit,{Dir => -1, Debug => 1});
      299 #
      300 #    print "testin $testin->{Name} $t_in  ---> testout $testout   \n";
      301 #
      302 #    my $testout2 = &Do_Transform($t_in,$bestfit,{Dir =>  1, Debug => 1});
      303 #
      304 #   print "testin $testin->{Name} $t_in  ---> testout2 $testout2   \n";
      305
      306 #    $out = &Do_Transform($im,$bestfit, {Dir => -1, Debug => 1});
      307 #    $out ->sethdr($h);
      308 #    wfits $out,'test1.fits';
      309 #
      310
      311
      312 #    my $test_in = pdl (10,500);
      313     my $test_in = pdl (910, 430);
      314     my $testout
      315         = &TrIm::Do_Transform(
      316             $test_in,
      317             $bestfit,
      318             {
      319                 Min_pixcoord_x => $min_x,
      320                 Max_pixcoord_x => $max_x,
      321                 Min_pixcoord_y => $min_y,
      322                 Max_pixcoord_y => $max_y,
      323                 Dir =>  -1, Debug => 1
      324             }
      325         );
      326     print "testin $test_in  ---> testout $testout   \n";
      327
      328     $out
      329         = &TrIm::Do_Transform(
      330             $out,
      331             $bestfit,
      332             {
      333                 Min_pixcoord_x => $min_x,
      334                 Max_pixcoord_x => $max_x,
      335                 Min_pixcoord_y => $min_y,
      336                 Max_pixcoord_y => $max_y,
      337                 Dir => +1, Debug => 1
      338             }
      339         );
      340     $out = float($out);
      341     $out ->sethdr($h);
      342
      343     return $out;
      344 }

    my $nord_x = 2;
    my $nord_y = 1;

    my $Nx = 1800;
    my $Ny = 1300;

    my $RA_0 = '17:13:44.282'; my $DEC_0 = '-37:06:15.09';