% ASG2Y_DEMO   A simple demonstration how to use
%    asg functions to simulate radiative transfer
%    for generated atmospheric states (by asg_demo)
%
%    Gives an example on basic usage of ASG2Y
%
%    The example simulate radiative transfer 
%    for a cloudy state generated by asg_demo 
%
%    A lot of default settings are made in this
%    demo, which may be re-considered depending
%    on what data should be used for
%
% FORMAT [G,y,ydata,dy]=asg2y_demo( Q, [ instrument ])
%
% OUT     G      structure with atmospheric variables on gformat data
%         y      As returned by *arts_y*.
%         ydata  As returned by *arts_y*.
%         dy     As returned by *arts_y*.        
%
% IN      Q       qarts setting structure
%
% OPT     sensor  choice of instrument 
%                 can be odin (default) or amsu      

% Created by Bengt Rydberg


function [G,y,ydata,dy]=asg2y_demo(Q,varargin)
%
[sensor] = optargs( varargin, { 'odin' } );

error( 'Can not be used presently. Not yet updated.' );


%-----------------------------------
%generate cloudy atmospheric states

[G,FS]=asg_demo( 3,sensor );
G=G{1};

%-----------------------------------
%SET Q-STUCTURE (for ARTS RT-simulations)

if ~exist('Q')
  %in sub-function set_q a lot of default settings are made
  Q=asg_q(sensor);
end
 

%------------------------------------
%ADD SSP OF ICE/WATER PARTICLES FOR THE RELEVANT INSTRUMENT FREQUENCIES 
FS.PROPS=asg_ssp(Q.F_GRID);
    
ind=find(strncmp({G.NAME},'Particles',9));
for j=ind
 G(j).PROPS=FS.PROPS(1).SSP(j-ind(1)+1);
end
  
%---------------------------------------
%RUN RTE CALCULATIONS
 
Q.P_GRID=G(1).GRID1;
Q.LAT_GRID=G(1).GRID2;
Q.LON_GRID=G(1).GRID3;

%adjust sensor_pos
if strcmp('amsu',sensor)
  
   %the sensor is placed right above the cloud 
   Q.SENSOR_POS(:,2)=FS.USE_LATS(1);
  
   %randomize a surface emissivity
   emis=min(0.95+0.05*randn(length(Q.F_GRID),1),1);
   emis=max(0.8,emis);
   FS.SURF_EMIS=emis;

end
  
if strcmp('odin',sensor)

   %make sure the sensor looks into the cloud field 
   if FS.LEG==1
      Q.SENSOR_POS(:,2)=-Q.SENSOR_POS(:,2)+FS.USE_LATS(1);
      Q.SENSOR_LOS(:,2)=180; 
   elseif FS.LEG==2
      Q.SENSOR_POS(:,2)=Q.SENSOR_POS(:,2)+FS.USE_LATS(1); 
      Q.SENSOR_LOS(:,2)=0; 
   else
      Q.SENSOR_POS(:,2)=-Q.SENSOR_POS(:,2)+FS.USE_LATS(1);
      Q.SENSOR_LOS(:,2)=180;
   end
  
end


%run arts simulations

workfolder =create_tmpfolder;

%don't consider water particles here/ should be updated
ind=find(strncmp({G.NAME},'LWC',3));
G=G(1:ind-1);

[y,dy,ydata]=asg_gfs2i(G,Q,FS,sensor,workfolder)

delete_tmpfolder(workfolder)





%-sub-functions


function Q=asg_q(instrument)

 Q = qarts;
 Q.OUTPUT_FILE_FORMAT='binary';

 %- Basic atmospheric variables
 %
 Q.ATMOSPHERE_DIM = 3;
 Q.P_GRID         = 100/16e3; % Unit here is pressure decade
 Q.LAT_GRID=[];
 Q.LON_GRID=[];
 Q.R_GEOID=6366700;

 %- Surface properties
 %
 if strcmp('amsu',instrument)
  se=0.6;
  Q.SURFACE_PROP_AGENDA = ...
  {'InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim, p_grid,lat_grid, lon_grid,rte_gp_p, rte_gp_lat, rte_gp_lon,t_field )', ...
  'Ignore(rte_los)',...
  ['NumericSet(surface_emissivity,',num2str(se),')'],'surfaceSimple'};
 end

 if strcmp('odin',instrument) 

   Q.SURFACE_PROP_AGENDA{1}=...
      ['InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim,'...
       'p_grid, lat_grid, lon_grid,rte_gp_p, rte_gp_lat,'...
       'rte_gp_lon,t_field)'];
   Q.SURFACE_PROP_AGENDA{2}='Ignore(rte_pos)';
   Q.SURFACE_PROP_AGENDA{3}='Ignore(rte_los)';
   Q.SURFACE_PROP_AGENDA{4}=['surfaceBlackbody(surface_los,'...
      'surface_rmatrix,surface_emission, f_grid, stokes_dim,surface_skin_t )'];
 
 end

 %- Spectroscopic stuff
 %
 Q.STOKES_DIM       = 1;
 Q.ABS_LINES_FORMAT = 'Arts';
 Q.ABSORPTION='OnTheFly';
 Q.ABS_LINESHAPE='Voigt_Kuntz6';
 Q.ABS_LINESHAPE_CUTOFF=-1;
 Q.ABS_LINESHAPE_FACTOR='quadratic';


 if strcmp('amsu',instrument)
  f_grid_file=fullfile( atmlab('ARTS_INCLUDES'),'amsu/amsub.f_grid_fast.xml');
  f_grid=xmlLoad(f_grid_file);
  Q.F_GRID           = f_grid;
  Q.ABS_LINES=fullfile(atmlab('ARTS_INCLUDES'),'amsu/amsub.hitran04_lines.xml');
  Q.ABS_MODELS={fullfile(atmlab('ARTS_INCLUDES'),'continua.arts')};
 end

 if strcmp('odin',instrument)
  Q.F_GRID           = [501.2e9 544.4e9];
  Q.ABS_LINES{1}=...
     fullfile(atmlab('ARTS_INCLUDES'),'odin/linefile.SM_AC2ab.xml'); 
  Q.ABS_LINES{2}=...
     fullfile(atmlab('ARTS_INCLUDES'),'odin/linefile.SM_AC1e.xml');
  Q.ABS_MODELS={fullfile(atmlab('ARTS_INCLUDES'),'odin/odinsmr.arts')};
 end

 %- RTE variables
 %
 Q.PPATH_STEP_AGENDA = { 'ppath_stepGeometric( ppath_step, atmosphere_dim, p_grid, lat_grid, lon_grid, z_field, r_geoid, z_surface, ppath_lmax)' };
Q.PPATH_LMAX=1e3;
Q.RTE_AGENDA= {'RteStd(iy,diy_dvmr,diy_dt,ppath,ppath_array,ppath_array_index,f_grid,stokes_dim,emission_agenda,abs_scalar_gas_agenda,rte_do_vmr_jacs,rte_do_t_jacs )'};
 Q.EMISSION_AGENDA={'emissionPlanck'};
 Q.IY_SPACE_AGENDA={'Ignore(rte_los)','Ignore(rte_pos)',...
              'MatrixCBR(iy,stokes_dim,f_grid)' };
 Q.J_DO=0;

 %- SENSOR variabels
 Q.ANTENNA_DIM=1;
 Q.Y_UNIT            = 'RJBT';
 Q.SENSOR_DO=0;


 %- Sensor POS and LOS
 %
 % Here hard coded
 %
 if strcmp('amsu',instrument)
   re=Q.R_GEOID;
   zplat = 850e3;
   xoff=[-7.5e3:2.5e3:7.5e3];
   yoff=[-7.5e3:2.5e3:7.5e3];
   xoff=repmat(xoff,length(xoff),1);
   yoff=repmat(yoff,1,length(yoff));
   xoff=xoff(:);
   yoff=yoff(:);
   s=sqrt(xoff.^2+yoff.^2);
   za=180-180/pi*atan(s/zplat);
   th=180/pi*atan(yoff./xoff);
   th(find(isnan(th)))=0;
   ivec=find(xoff==0 & yoff<0);
   th(ivec)=th(ivec)+360;
   ivec=find(yoff==0 & xoff<0);
   th(ivec)=180;
   ivec=find(xoff<0 & yoff>0);
   th(ivec)=th(ivec)+180;
   ivec=find(xoff<0 & yoff<0);
   th(ivec)=th(ivec)+180;
   ivec=find(xoff>0 & yoff<0);
   th(ivec)=th(ivec)+360;
   ivec=find(th>180);
   th(ivec)=th(ivec)-360;
   ppos=zeros(length(th),3);
   ppos(:,1)=re+zplat;
   plos=zeros(length(th),2);
   plos(:,1)=za;
   plos(:,2)=th;
   Q.SENSOR_POS=ppos;
   Q.SENSOR_LOS=plos;
 end

 if strcmp('odin',instrument)
   re=Q.R_GEOID;
   zplat = 600e3;
   tan=[-5e3 -2e3:5e2:-500 0:2.5e2:10e3 10.5e3:5e2:12e3 15e3]';
   tan=6e3;
   za = geomztan2za(re,zplat,tan);
   ppos=zeros(length(za),3);
   ppos(:,1)=re+zplat;
   ppos(:,2)=-21.5;
   plos=zeros(length(za),2);
   plos(:,1)=za;
   Q.SENSOR_POS=ppos;
   Q.SENSOR_LOS=plos;
 end

 %- MC variables
 %
 Q.CLOUDBOX.METHOD                      = 'MC';
 Q.CLOUDBOX.METHOD_PRMTRS.STD_ERR       = 3;      
 Q.CLOUDBOX.METHOD_PRMTRS.MAX_TIME      = -1;
 Q.CLOUDBOX.METHOD_PRMTRS.MAX_ITER      = -1;
 Q.CLOUDBOX.METHOD_PRMTRS.Z_FIELD_IS_1D = 0;
 Q.CLOUDBOX.LIMITS                      = [];
 Q.CLOUDBOX.OPT_PROP_GAS_AGENDA={'ext_matInit','abs_vecInit',...
                                'ext_matAddGas','abs_vecAddGas'};
 Q.CLOUDBOX.OPT_PROP_PART_AGENDA={'ext_matAddPart','abs_vecAddPart'};
 Q.CLOUDBOX.SPT_CALC_AGENDA={};
 Q.CLOUDBOX_DO=0;

return

%-------------------------------------------------------------------------

function P=asg_ssp(f_grid)
 %properties for adding pnd fields
 alpha=0;Ni=10;xnorm=50e-6;
 [x_i,w_i]=gauss_laguerre(alpha,Ni,xnorm);
 Nj=3;xnormw=10e-6;
 [x_k,w_k]=gauss_laguerre(alpha,Nj,xnormw);
 %x_i will be the diameter of the particles considered
 %calculate ssp
 for i=1:2

  if i==1 
   F=create_ssp('sphere',f_grid,x_i/2);
   P(i).x=x_i;
   P(i).w=w_i;
   P(i).x_norm=xnorm;
   P(i).PSD='MH97';
  else
   T_grid=[273 283 293 303 313];
   for i1=1:length(f_grid)
     for j1=1:length(T_grid)
       rfr_index(i1,j1) = sqrt(eps_water_liebe93( f_grid(i1), T_grid(j1) ));
     end
   end
   F=create_ssp('sphere',f_grid,x_k,T_grid,rfr_index);
   P(i).x=x_k;
   P(i).w=w_k;
   P(i).x_norm=xnormw;
   P(i).PSD='Water';
  end
  P(i).SSP=F;
  P(i).method='gauss-laguerre';
  P(i).shape='sphere';  
 end

return


%-------------------------------------------------------
function [y,dy,ydata]=asg_gfs2i(G,Q,FS,sensor,workfolder)

 %- Determine necessary size of cloud box and adjust it
 %
 [Q,G]=asg_set_cb(Q,G);


 %check if there is any liquid cloud
 %if not we don't need to include it in the simulations
 ind=find(strcmp({G.NAME},'LWC field'));
 if ~isempty(ind)
   ind=find(strcmp({G.NAME},'LWC field'));
   G(ind).PROPS=['liquidcloud-MPM93'];
   G(ind).DATA=G(ind).DATA*1e-3;
   G(ind).DATA_NAME='Volume mixing ratio';
   G(ind)= asg_regrid( asgD, G(ind), Q );
   G(ind).DATA(find(G(ind).DATA(:)>10e-3))=10e-3;
 else
   ind=setdiff(1:length(G),ind);
   G=G(ind);
 end

 %- Run ARTS
 %
 if strcmp(sensor,'amsu')
   

  %SURFACE EMISSIVITIES
  Se=[];
  for j=1:length(FS.SURF_EMIS)
    se=num2str(FS.SURF_EMIS(j));
    Se=[Se,',',se];
  end         
  Se=Se(2:end);

  Q.SURFACE_PROP_AGENDA= ...
    {['InterpAtmFieldToRteGps(surface_skin_t,atmosphere_dim,'...
      'p_grid, lat_grid, lon_grid,rte_gp_p, rte_gp_lat,'...
      'rte_gp_lon,t_field)'],...
      'Ignore(rte_pos)',...
      ['VectorCreate(surface_emissivity)'],... 
     ['VectorSet(surface_emissivity,[',Se,'])'],...
      ['surfaceFlatVaryingEmissivity(surface_los, surface_rmatrix,'...
      'surface_emission, f_grid, stokes_dim,atmosphere_dim, rte_los,'...
      'surface_skin_t,surface_emissivity)']};  
   Q = asg2q( asgD, G, Q, workfolder );
   Q.ABS_NLS=[];
   %create an include file
   fid=fopen(fullfile(workfolder,'include.arts'))
   if fid==-1
     fid=fopen(fullfile(workfolder,'include.arts'),'a');
     fprintf(fid,'Arts2 {\n'); 
     fprintf(fid,'#VectorSetExplicitly(abs_n2,[0.7808])\n');
     fprintf(fid,'VectorSet( abs_n2, [ 0.7808 ] )\n');
     fprintf(fid,'IndexSet( abs_p_interp_order,   5 )\n');
     fprintf(fid,'IndexSet( abs_t_interp_order,   7 )\n');
     fprintf(fid,'IndexSet( abs_nls_interp_order, 5 )}\n');
     fclose(fid);
   end
   Q.INCLUDES= {fullfile(workfolder,'include.arts')};
   Q                  = qarts_abstable( Q,9,100,[],[] );
   [y,ydata,dy] = arts_y( Q, workfolder );

 end


 if strcmp(sensor,'odin')

  for i= 1:length(Q.F_GRID)
    if 0
     if i==1
       cind=find(strncmp(lower({G.NAME}),'hno3',4));
       ind=setdiff(1:length(G),cind);
       G1=G(ind);
     else
       cind=find(strncmp(lower({G.NAME}),'clo',3));
       ind=setdiff(1:length(G),cind);
       G1=G(ind);
     end
   else
    G1=G;
   end
    
    Q1           = asg2q( asgD, G1, Q, workfolder );
    Q1.ABS_NLS    =[];
    Q1           = qarts_abstable( Q1,9,100,[],[] );
    
    %create an include file
    fid=fopen(fullfile(workfolder,'include.arts'))
    if fid==-1
     fid=fopen(fullfile(workfolder,'include.arts'),'a');
     fprintf(fid,'Arts2 {\n'); 
     fprintf(fid,'#VectorSetExplicitly(abs_n2,[0.7808])\n');
     fprintf(fid,'VectorSet( abs_n2, [ 0.7808 ] )\n');
     fprintf(fid,'IndexSet( abs_p_interp_order,   5 )\n');
     fprintf(fid,'IndexSet( abs_t_interp_order,   7 )\n');
     fprintf(fid,'IndexSet( abs_nls_interp_order, 5 )}\n');
     fclose(fid);
    end
    Q1.INCLUDES= {fullfile(workfolder,'include.arts')};
    Q1.ABS_LINES =Q.ABS_LINES{i};
    Q1.F_GRID    =Q.F_GRID(i);
    [y(:,i),ydata{i},dy(:,i)] = arts_y( Q1, workfolder );
  
  end

 end


return

%-----------------------------------------------------------------

function [Q,G]=asg_set_cb(Q,G)

iwc_ind = find( strcmp( lower({G.NAME}), 'iwc field' ) );
[C1,C2,C3,C4,C5,C6]=deal([]);

for i=1:size(G(iwc_ind).DATA,3)
  [c1,c2]=find(G(iwc_ind).DATA(:,:,i));
  if ~isempty(c1)
     C1=min([C1,vec2col(c1)']);
     C2=max([C2,vec2col(c1)']);
  end
  if ~isempty(c2)
     C3=min([C3,vec2col(c2)']);
     C4=max([C4,vec2col(c2)']);
  end
end


for i=1:size(G(iwc_ind).DATA,1)
  [c1,c2]=find(squeeze(G(iwc_ind).DATA(i,:,:)));
    if ~isempty(c2)
       C5=min([C5,vec2col(c2)']);
       C6=max([C6,vec2col(c2)']);
   end
end
%

if ~isempty(C1) %if C1 is empty we don't have any cloud particles
    p1   = G(iwc_ind).GRID1(C1);
    p2   = G(iwc_ind).GRID1(C2);
    lat1 = G(iwc_ind).GRID2(C3);
    lat2 = G(iwc_ind).GRID2(C4);
    z    = p2z_cira86([p1 p2]',mean([lat1 lat2]),160);
    lon1 = G(iwc_ind).GRID3(C5);
    lon2 = G(iwc_ind).GRID3(C6);
    %
    Q.CLOUDBOX.LIMITS = [z(1)-2e3 z(2)+2e3 lat1 lat2 lon1 lon2];
    
else
    %remove all particles in G 
    ind  =find(strncmp({G.NAME},'Particles',9)==0);
    G    = G(ind);

end