include 'image.g';

#plots a horizon scan assuming a given integration has no rfi (intref)
rfihorizonplot_ref := function(msname,scan,intref1,intref2,numif,numpol,imname,feed=1,azpixsize=1.5) {

    d.open(msname);

    r1 := d.getr(scan,1,pol=numpol,nif=numif,int=intref1,nfeed=feed);
    r2 := d.getr(scan,1,pol=numpol,nif=numif,int=intref2,nfeed=feed);

    # get range of integrations to do
    ints := d.qscan(scan).ints;

    data := array(0,1,len(r1.data.arr[1,]),ints,1,1);

    slope := r1;
    slope.data.arr := (r2.data.arr-r1.data.arr)/(intref2-intref1);
    intercept := r1;
    intercept.data.arr := r1.data.arr-slope.data.arr*intref1;

    for (i in 1:ints) {
	r := r1;
	r.data.arr := intercept.data.arr+slope.data.arr*i;
	s := d.getr(scan,1,pol=numpol,nif=numif,int=i,nfeed=feed);
	for (j in 1:len(s.data.arr[1,])) {
	    where := len(r.data.desc.chan_freq.value)+1-j;
	    data[1,j,i,1,1] := s.data.arr[1,where]-r.data.arr[1,where];
	}
    }


    #create coordinate system
    cs := coordsys(linear=1,spectral=T,stokes='I',direction=T);
    cs.reorder([2,3,4,1]);
    cs.setobserver('T. Minter');
    ep := dm.epoch('UTC', 'today')
    cs.setepoch(ep);
    cs.settelescope('GBT');
    freqs := r.data.desc.chan_freq;
    for (j in 1:len(r.data.desc.chan_freq.value)) {
	where := len(r.data.desc.chan_freq.value)+1-j;
	freqs.value[j] := r.data.desc.chan_freq.value[where];
    }
    cs.setspectral(frequencies=freqs);
    cs.setrestfrequency(r.data.desc.restfrequency);
    names := cs.names();
    names[3] := 'Azimuth';
    cs.setnames(names);
    increment := cs.increment();
    increment[3]:=azpixsize;
    cs.setincrement(increment);
    cs.setunits("deg","linear",T);
    cs.setreferencecode('TOPO','spectral',T);
    cs.summary();

    #make image and display
    myim := imagefromarray(imname,data,cs,T,T);
#    myim.view();

    print imname;
#    bogus := readline();
    myim.done();

    d.close(msname);

}

#plots a elevation scan assuming a given integration has no rfi (intref)
rfitipplot_ref := function(msname,scan,intref1,intref2,numif,numpol,imname,feed=1,
		       elpixsize=0.5833333333333333333333333,way=1) {

    d.open(msname);

    r1 := d.getr(scan,1,pol=numpol,nif=numif,int=intref1,nfeed=feed);
    r2 := d.getr(scan,1,pol=numpol,nif=numif,int=intref2,nfeed=feed);

    # get range of integrations to do
    ints := d.qscan(scan).ints;

    data := array(0,1,len(r1.data.arr[1,]),ints,1,1);

    slope := r1;
    slope.data.arr := (r2.data.arr-r1.data.arr)/(intref2-intref1);
    intercept := r1;
    intercept.data.arr := r1.data.arr-slope.data.arr*intref1;

    for (i in 1:ints) {
	r := r1;
	r.data.arr := intercept.data.arr+slope.data.arr*i;
	s := d.getr(scan,1,pol=numpol,nif=numif,int=i,nfeed=feed);
	for (j in 1:len(s.data.arr[1,])) {
	    where := len(r.data.desc.chan_freq.value)+1-j;
	    data[1,j,i,1,1] := s.data.arr[1,where]-r.data.arr[1,where];
	}
    }


    #create coordinate system
    cs := coordsys(linear=1,spectral=T,stokes='I',direction=T);
    cs.reorder([2,3,4,1]);
    cs.setobserver('T. Minter');
    ep := dm.epoch('UTC', 'today')
    cs.setepoch(ep);
    cs.settelescope('GBT');
    freqs := r.data.desc.chan_freq;
    for (j in 1:len(r.data.desc.chan_freq.value)) {
	where := len(r.data.desc.chan_freq.value)+1-j;
	freqs.value[j] := r.data.desc.chan_freq.value[where];
    }
    cs.setspectral(frequencies=freqs);
    cs.setrestfrequency(r.data.desc.restfrequency);
    names := cs.names();
    names[3] := 'Elevation';
    cs.setnames(names);
    if (way != 1) {
	step := -1*elpixsize;
	rv := 80.;
    } else {
	step := elpixsize;
	rv := 10.;
    }
    cs.setreferencevalue(type='linear',value=rv);
    increment := cs.increment();
    increment[3]:=step;
    cs.setincrement(increment);
    cs.setunits("deg","linear",T);
    cs.setreferencecode('TOPO','spectral',T);
    cs.summary();

    #make image and display
    myim := imagefromarray(imname,data,cs,T,T);
#    myim.view();

    print imname;
#    bogus := readline();
    myim.done();

    d.close(msname);

}


#plots a horizon scan assuming rfi is present in all scans
rfihorizonplot := function(msname,scan,numif,numpol,imname,feed=1,
			   azpixsize=1.5,removebaseline=T,nsize=25) {

    d.open(msname);

    # get range of integrations to do
    ints := d.qscan(scan).ints;

    r := d.getr(scan,1,pol=numpol,nif=numif,int=intref1,nfeed=feed);
    data := array(0,1,len(r.data.arr[1,]),ints,1,1);


    

    for (i in 1:ints) {

	if (i%10==0) {
	    printf('%d/%d...',i,ints);
	}

	ss := d.getr(scan,1,pol=numpol,nif=numif,int=i,nfeed=feed);
	s := ss;
	if (removebaseline) {
	    s.data.arr[1,] -:= rfibaseline(ss,nsize).data.arr[1,];
	}
	for (j in 1:len(s.data.arr[1,])) {
	    data[1,j,i,1,1] := s.data.arr[1,j];
	}
    }

    #create coordinate system
    cs := coordsys(linear=1,spectral=T,stokes='I',direction=T);
    cs.reorder([2,3,4,1]);
    cs.setobserver('T. Minter');
    ep := dm.epoch('UTC', 'today')
    cs.setepoch(ep);
    cs.settelescope('GBT');
    freqs := r.data.desc.chan_freq;
    for (j in 1:len(r.data.desc.chan_freq.value)) {
	where := len(r.data.desc.chan_freq.value)+1-j;
	freqs.value[j] := r.data.desc.chan_freq.value[where];
    }
    cs.setspectral(frequencies=freqs);
    cs.setrestfrequency(r.data.desc.restfrequency);
    names := cs.names();
    names[3] := 'Azimuth';
    cs.setnames(names);
    increment := cs.increment();
    increment[3]:=azpixsize;
    cs.setincrement(increment);
    cs.setunits("deg","linear",T);
    cs.setreferencecode('TOPO','spectral',T);
    cs.summary();

    #make image and display
    myim := imagefromarray(imname,data,cs,T,T);
#    myim.view();

    print imname;
#    bogus := readline();
    myim.done();

    d.close(msname);

}


# take the tenth lowest values of the data within nsize data points in each
# direction as the baseline value
rfibaseline := function(data,nsize=25) {

    if (nsize <10) {
	print "nsize should be greater than 9";
	return;
    }

    newdata := data;

    maxsize := len(data.data.arr[1,]);

    for (i in 1:maxsize) {

	lo_end := i-nsize;
	hi_end := i+nsize;
	if (lo_end < 1) lo_end:=1;
	if (hi_end > maxsize) hi_end:=maxsize;

	mydata := data.data.arr[1,lo_end:hi_end];
	mydata := sort(mydata);
	newdata.data.arr[1,i] := mydata[10]; 
    }

    return newdata;
}