;  POLAR_MEMDEMO -- a brief script to demonstrate operation of
;    hsi_polar_mem, which takes a map, converts it to fake observed 
;    visibilities
;    and then shifts them to Sun centre (lifted from EJS's polar_mapper
;    code (Nov. 23, version 4). A MEM reconstruction is then performed
;    to reconstruct a map from these visibilities.
;
;  AIM -- To reconstruct a map from simulated observations of a specified model
;    map using MEM in polar co-ords.
;
;  VERSION -- This is designed to work with hsi_polar_mem Vers. 4
;  RESULTS -- TBA
;
;  CONCLUSION -- TBA
;
;  TO DO -- agreement of estimate vs actual visi errors is very good
;           if xypos=[0,0] for a single gaussian point source.
;           For other xypos values, the result is variable:
;           some are good, some involve the actual~ 2*estimated.
;           I believe this is to do with polar pixel boundaries
;           vs. the pixels used in the simulation. It could also be
;           something to do with binning, in which case a revision
;           of pmap2vis's binning procedures are needed.
;           From RAS: Gaussians, not pixels, are used in the simulation.
;           note above could still be relevant with errors creeping in
;           because of "pixelisation" of assumed pointed sources.  
;           
;           Significant disgreement still exists for det=3 (and possibly <3)
;            my guess was this is because the number of points per cycle
;            is getting quite low. The fact that the disagreement
;            seems worst near the limb is consistent with this.
;            HOWEVER, inspecting plots shows that there is 
;            a noticeable phase difference between noisy and perfect visis.
;            EVEN when xypos=0.0.
;            Perhaps it is most easily seen in the errors for the finer dets 
;            because they are naturally most sensitive to phase shifts.
;            When xypos are nonzero, the errors are definitely because
;            if phase shifts. It is possible that all the above problems
;            have the same underlying cause: pixel boundary and binning
;            causing systematic shifts in phase.
;
;	     
;           There is disagreement for det=8
;            I believe this is due to the very long cycle involved. 
;
;    AJC -- Feb 22, 2000
;


  MSflag=0 ; set this to 1 to run mem sato, -1 to set up only

  sim_photons=10000 ; Number of photons/s/SC
  heap_gc
  o=obj_new('hsi_image') 
       
  timebindef = o->get( /time_bin_def ) 
  timebinmin=512

  btot=4l*(sim_photons/(2l^22/(timebindef[8]*timebinmin)))

  ng=1
  gaussians = {gaussian_source_str}
  gaussians = replicate( {gaussian_source_str}, ng )
      
  gaussians[0].amplitude =1.0
  gaussians[0].xypos(0)  = 0.5
  gaussians[0].xypos(1)  = 0.0
  gaussians[0].xysigma  = [0.1,0.1]       
;  gaussians[1].amplitude =1.0
;  gaussians[1].xypos(0)  = 40.0
;  gaussians[1].xypos(1)  = 0.0
;  gaussians[1].xysigma  = [5,10]    
;  gaussians.tilt_angle_deg= 0.0

       RHO_MAP=600.0 & AZ_MAP=0.043945312 ; location of map center in polars
       xy_offset=[RHO_MAP*cos(-AZ_MAP*!dtor),RHO_MAP*sin(-AZ_MAP*!dtor)];in xy

       sim_xy_offset=xy_offset+0 ; Is this an error in the flare posn? ASK EJS
       det_mask=intarr(27) 
;**********
       for det_index=0,8 do det_mask(det_index)=1 
;       det_mask[3]=1
;**********			     
       time_range=[0,4]  ; If this is changed, polar_mapper won't work
       pixel_size=[1.0,1.0]
       sim_pixel_size=pixel_size[0]
       sim_model=gaussians
       M_pmap=hsi_gauss2pmap(sim_model,btot,sim_pixel_size,xy_offset,/verbose)
  
       if(MSflag eq 1) then begin
         t1=systime(1)
         o->set,$
             SIM_TIME_range=time_range, sim_a2d=det_mask,$
             TIME_RANGE=time_range, SIM_PHOT=sim_photons, $
             DET_INDEX_MASK=det_mask,TIME_BIN_DEF=timebindef, $
             TIME_BIN_MIN=timebinmin,SIM_MODEL=sim_model, $
             SIM_XYOFFSET=sim_xy_offset, XYOFFSET=xy_offset,$   
             PIXEL_SIZE=pixel_size,SIM_PIXEL_SIZE=sim_pixel_size,image_algorithm='mem'
         map=o->getdata()
         print,'Time taken for MEM= ',systime(1)-t1
       endif else if msflag eq 0 then begin
         t1=systime(1)
         o ->set,$
             SIM_TIME_range=time_range, sim_a2d=det_mask,$
             TIME_RANGE=time_range, SIM_PHOT=sim_photons, $
             DET_INDEX_MASK=det_mask,TIME_BIN_DEF=timebindef, $
             TIME_BIN_MIN=timebinmin,SIM_MODEL=sim_model, $
             SIM_XYOFFSET=sim_xy_offset, XYOFFSET=xy_offset,$   
             PIXEL_SIZE=pixel_size,SIM_PIXEL_SIZE=sim_pixel_size,image_algorithm='polar'
         map=o->getdata()
         print,'Time taken for polar MEM= ',systime(1)-t1
       endif
end