guisdap Code

function [varargout]=plasma_summary(pl_dir,printdir,expt,gates,plots)
% function plasma_summary(pl_dir,printdir,expt,gates,plots)
% GUISDAP 8.6  8 Oct 2010  Copyright EISCAT
% Get/plot summary for plasmalines
% Using setup in pl_def.m in experiment directory
% input: pl_dir Directory containing integrated plasma lines (needed)
%        printdir Directory to save plot to (eps and png)
%        expt   Name of experiment (or guessed from pl_dir)
%        gates  Gates to display
%        plots  Display a number of plots for internal checks (-1, don't clear)
% output: pl structure containing
%            .t Time
%            .s Spectra
%            .f Freqeuencies
if nargin<2, printdir=[]; end
if nargin<3, expt=[]; end
verbose=0;
if nargin<4
 gates=[];
elseif gates==0
 gates=[]; verbose=1;
end
gate=[];
if nargin<5, plots=0; end
%if isempty(plots), plots=0; end
if plots~=-1, clf(findobj('type','figure','userdata',7)), plots=0; end
if isempty(expt)
 [i,expt]=fileparts(pl_dir);
 if isempty(expt), [i,expt]=fileparts(i); end
end
global path_exps path_GUP local
%Get plasmaline data
[fl,msg]=getfilelist(fullfile(pl_dir,filesep));
if msg, error(msg), end
nfl=length(fl);
load(canon(fullfile(fl(1).dir,sprintf('%08d%s',fl(1).file,fl(1).ext))))

exps=dir(path_exps); pldef='pl_def'; pdefs=[];
for i=1:length(exps)
 d=fullfile(path_exps,exps(i).name,pldef);
 if exist([d '.m'],'file')
  pdefs=[pdefs i];
  if strfind(expt,exps(i).name) & exist([d '.m'],'file')
   run(d),break
  end
 end
end

if ~exist('startad','var')
 pdefs=sprintf('%s ',exps(pdefs).name);
 warning(['No such experiment defined ( ' pdefs ')'])
 verbose=1;
 ran=[200 300]; nfft=128; nint=1; ngates=1; nlag=0; maxe=2; nup_d=1; skip_if=0;
 freq=[-5 5]*1e6; dt=1e-6; invert=1; fradar=930e6; ele=77.5; updown=0:1;
 startad=0;
end
if isempty(gates), gates=1:ngates; end
if verbose
 ran=minput('Ranges (km)',ran);
 nfft=minput('No fft points',nfft);
 nint=minput('No to integrate',nint);
 ngates=minput('No of gates',ngates);
 nlag=minput('No of lags',nlag);
 maxe=minput('Max sigma',maxe);
 nup_d=minput('No frequencies',nup_d);
 skip_if=minput('Skip interference',skip_if);
 freq=minput('Offset frequencies',freq);
 dt=minput('Sampling frequency',dt);
 invert=minput('Invert frequency scale',invert);
 fradar=minput('Radar frequency',fradar);
 ele=minput('Approx elevation',ele);
 updown=minput('Up/Down shifted (0/1)',updown);
 startad=minput('Start address',startad);
 gates=minput('Gates to display',gates);
 upar0=minput('Start of upar',46);
else
 upar0=46;
end
nfreq=length(freq); dgates=length(gates);
if ~ismember(gate,gates), gate=gates(end); end

endad=startad+(nfft+nlag)*ngates*nint*nfreq/length(startad)-1;
add=[]; for i=1:length(startad), add=[add startad(i):endad(i)]; end
if nlag>0, nfft=2*nlag; end
plspec=zeros(dgates,nfft,nfreq,nfl);
p_time=zeros(2,nfl);
upar=zeros(nfreq,nfl);
for i=1:nfl
 if i>1
  filename=fullfile(fl(i).dir,sprintf('%08d%s',fl(i).file,fl(i).ext));
  load(canon(filename,0))
 end
 p_time(:,i)=datenum(row(d_parbl(1:6)))+[-d_parbl(7)/86400;0];
 if nlag>0
  d=sum(reshape(d_data(add),ngates,nlag,nint,nfreq),3)/d_parbl(22);
  d=reshape(d(gates,:,1,:),dgates,nlag,nfreq);
  d=real(fft([d zeros(dgates,1,nfreq) conj(d(:,nlag:-1:2,:))],[],2));
  d=d(:,[nfft/2+1:nfft 1:nfft/2],:);
 else
  d=sum(reshape(real(d_data(add)),nfft,ngates,nint,nfreq),3)/d_parbl(22);
  d=reshape(d([nfft/2+1:nfft 1:nfft/2],gates,1,:),nfft,dgates,nfreq);
  d=permute(d,[2 1 3]);
 end
 plspec(:,:,:,i)=d;
 if exist('uparfreq','var')
  upar(:,i)=d_parbl(upar0+(1:nfreq));
 end
end
if plots
 imagesc(reshape(plspec(gate,:,:,:),nfft*nfreq,nfl)), pause
end
plsig=zeros(size(plspec));

bands=unique(upar(1,:));
for band=1:length(bands)
 b=find(upar(1,:)==bands(band)); lb=length(b);
 %Get spectral shape and find interferences
 disp('Removing noise')
 wrap=9;
 for g=1:dgates
  spec_shape=reshape(median(plspec(g,:,:,b),4),nfft,nfreq);
  filt_shape=spec_shape;
  fsh=[filt_shape(end-wrap+1:end,:);filt_shape;filt_shape(1:wrap,:)];
  d2=diff(filt_shape);
   fd2=2*std(d2); jj=-(wrap-1):nfft+wrap;
  for i=1:nfreq
   dd=find(d2(:,i)<-fd2(i)); du=find(d2(:,i)>fd2(i));
   d=unique([dd;du+1]); jd=jj; jd(d+wrap)=[];
   filt_shape(:,i)=interp1(jd,fsh(jd+wrap,i),1:nfft,'spline');
  end
  if plots
   plot([filt_shape spec_shape]), pause
  end
  interf_spec=spec_shape-filt_shape;
  if skip_if
   interf_spec(find(interf_spec~=0))=NaN;
  %filt_shape=spec_shape;
  end
 
 %Find a timevarying background (Use Tsys instead?)
  fp=reshape(median(plspec(g,:,:,b),2),nfreq,lb);
  nplspec=zeros(nfft,nfreq,lb);
  spnorm=median(fp,2);
  %spnorm=median(filt_shape,1);
  for i=1:nfreq
   nplspec(:,i,:)=filt_shape(:,i)*fp(i,:)/spnorm(i);
  end
  %Remove background and interference
  plsig(g,:,:,b)=reshape(plspec(g,:,:,b),nfft,nfreq,lb)-nplspec-reshape(interf_spec(:)*ones(1,lb),nfft,nfreq,lb);
  %Normalise signal
  %plsig(g,:,:,b)=plsig(g,:,:,b)./reshape(nplspec,1,nfft,nfreq,lb);
 end
end
if plots
 imagesc(reshape(plsig(gate,:,:,:),nfft*nfreq,[])), pause
end
 
freq_scale=invert*(-nfft/2:nfft/2-1)/dt/nfft;
if exist('uparfreq','var')
 if isnan(uparfreq(1))
  freq=1e6*upar;
 else
  freq=(uparfreq'*ones(1,nfl)+invert*upar)*1e6;
 end
else
 freq=col(freq)*ones(1,nfl);
end
df=diff(freq_scale(1:2));
plfreq=reshape(freq_scale'*ones(1,nfreq*nfl)+ones(nfft,1)*reshape(freq,1,nfl*nfreq),nfft,nfreq,nfl);
if nargout
 varargout(1)={struct('t',p_time,'f',plfreq,'s',plsig)};
 if nargout>1
  varargout(2)={struct('ran',ran,'ele',ele,'updown',updown,'nup_d',nup_d,'df',df,'fradar',fradar,'maxe',maxe,'gate',gate)};
 end
end
%plsig=plsig/median(plsig(:));
%plsig(find(plsig<1))=1; plsig=log(plsig);
[d,forder]=sort(freq,1);
for i=1:nfl
 plsig(:,:,:,i)=plsig(:,:,forder(:,i),i);
 plfreq(:,:,i)=plfreq(:,forder(:,i),i);
end
if plots
 plot(plfreq(:,:,round(nfl/2)),reshape(plsig(gate,:,:,round(nfl/2)),nfft,nfreq))
 pause
elseif plots==0
 npanel=dgates*nfreq;
 height(2)=.03; height(1)=(0.85-(npanel-1)*height(2))/npanel;
 set(gcf,'Position',[200 30 587 807],'renderer','painters','PaperPosition',[0.4 0.7 20.65 28.4],'UserData',7)

 for g=1:dgates
  mmmsig=mean(max(max(plsig(g,:,:,:))));
  for j=1:nfreq
   pn=((g-1)*nfreq+j-1);
   plfreqj=reshape(([plfreq(:,j,:);plfreq(end,j,:)+df]-df/2)/1e6,nfft+1,nfl);
   axes('Position',[.12 .06+pn*sum(height) .85 height(1)])
   for i=1:nfl
    surface(p_time(:,i),plfreqj(:,i),[plsig(g,:,j,i) plsig(g,end,j,i)]'*ones(1,2))
   end
   shading flat, caxis([0 mmmsig])
   ymax=max(max(plfreqj));
   if pn==0, xlabel('Time [UT]'), end
   mydatetick(gca,p_time([1 end]),pn==0)
   set(gca,'ylim',[min(min(plfreqj)) ymax],'box','on','tickdir','out','xgrid','on','ygrid','on')
   if j==1
    text('Position',[p_time(end) ymax],'VerticalAlignment','bottom','FontSize',12,...
    'HorizontalAlignment','Right','String',sprintf('Range %.0f-%.0f km',ran(gates(g),:)))
   else
    set(gca,'xticklabel',[])
   end
  end
 end
 [h,d]=strtok(d_ExpInfo);
 colormap(vizu('myb'))
 %load(fullfile(path_GUP,'matfiles','logo'))
 %axes('Position',[.07 .92 .08 .06]); plot(y,x,'.k')
 %set(get(gca,'child'),'markersize',1)
 %set(gca,'xlim',[0 202],'ylim',[0 202],'visible','off')
 axes('Position',[.07 .92 .08 .06]); eiscatlogo(.5)
 text('Position',[11 10],'VerticalAlignment','top','FontSize',16,...
  'HorizontalAlignment','Left','FontWeight','bold',...
  'String','EISCAT Scientific Association','UserData','Copyright');
 text('Position',[11 -8],'interpreter','none',...
  'String',sprintf('Plasma lines  %s  %s',d,datestr(mean(p_time(1,:)),1)),'UserData','Experiment')
 text('Position',-[15 150],'String','Frequency offset [MHz]','Rotation',90,...
  'VerticalAlignment','middle','HorizontalAlignment','center')
 if ~isempty(printdir)
  if pl_dir(end)==filesep, pl_dir(end)=[]; end
  if pl_dir=='.', pl_dir=pwd; end
  [dum,fname]=fileparts(pl_dir);
  fname=minput('Print file name (.pdf .png)',fullfile(printdir,[fname '_plasmaline']),1);
  if isunix & local.x
    gd=fullfile(matlabroot,'sys','ghostscript',filesep);
    gsbin=fullfile(gd,'bin',lower(computer),'gs');
    gsinc=sprintf('-I%sps_files -I%sfonts',gd,gd);
    if ~exist(gsbin,'file'), gsbin='gs'; gsinc=[]; end
    print(gcf,'-opengl','-depsc2','-r600',[fname '.eps'])
    unix(sprintf('%s %s -dNOPAUSE -dFitPage -q -sDEVICE=pdfwrite -sOutputFile=%s.pdf %s.eps </dev/null >/dev/null',gsbin,gsinc,fname,fname));
    unix(sprintf('%s %s -dNOPAUSE -dFitPage -q -sDEVICE=png256 -sOutputFile=%s.png %s.eps </dev/null >/dev/null',gsbin,gsinc,fname,fname));
    delete([fname '.eps'])
   %print(gcf,'-opengl','-dpdf','-r600',[fname '.pdf'])
   %print(gcf,'-dpng256',[fname '.png'])
  else
   print(gcf,'-dpdf',[fname '.pdf'])
   print(gcf,'-dpng',[fname '.png'])
  end
  fprintf('%s.pdf and .png saved\n',fname);
  insert_exif(gcf,fname,{'pdf' 'png'})
 end
end