% Tsys, TxPw  on 991025
%******************************* OBS !! *******************************
% This is presented by Saito_san & Ogawa_san !! and arranged by Sugino.
% Please use ionomodel.m which declears global variables 'MCA' and 'MRA',
% and save_results.m which saves variables 'MCA' and 'MRA'
% in analyzing by guisdap.
% gc2gg.m and gg2gc.m enable to calculate latitude and longtitude.
% These originals are in /work11/saito11/m165_dir.
% The direction of Vi analyzed by GUISDAP is opposite 
% to that by EISCAT standard which is positive for going away.
% You can choose both expression by rewriting make_dasx.m files.
% Default is positive for coming closer.
%********************************************************************** 


%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& Edit here &&&&&&&&&&&&&&&&&&&&&&&&&&&&&
% date
sdate0='920330';
name_expr_set='CP1K'
name_site='T';
result_path='/EISCAT00/920330/20sec/tromso/'
% directory of output file (.gdas)
    output_path=  result_path
% direction of Vi,   GUISDAP default:0,   EISCAT standard:1
    veloc=1;
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

if name_site == 'T'
 site0='Tromso';
 ss='tr';
 site_lat=69.+35/60.+11.239/3600.;   % r_RECloc(1)
 site_lon=19.+13/60.+38.746/3600.;   % r_RECloc(2)
 site_lev=0.0865;   % r_RECloc(3)
end
if name_site == 'K'
 site0='Kiruna';
 ss='ki';
 site_lat=67.+51/60.+38.373/3600.;   % r_RECloc(1)
 site_lon=20.+26/60.+06.806/3600.;   % r_RECloc(2)
 site_lev=0.4179;   % r_RECloc(3)
end
if name_site == 'S'
 site0='Sodankyla';
 ss='so';
 site_lat=67.+21/60.+49.266/3600.;   % r_RECloc(1)
 site_lon=26.+37/60.+37.497/3600.;   % r_RECloc(2)
 site_lev=0.1973;   % r_RECloc(3)
end
if name_site == 'V'
 site0='Tromso VHF';
 ss='tr';
 site_lat=69.+35/60.+11.9408/3600.;   % r_RECloc(1)
 site_lon=19.+13/60.+13.2300/3600.;   % r_RECloc(2)
 site_lev=0.0853;   % r_RECloc(3)
end
if name_site == 'L'
 site0='Longyearbyen';
 ss='lo';
 site_lat=78.+09/60.;   % r_RECloc(1)
 site_lon=16.+02/60.;   % r_RECloc(2)
 site_lev=0.434;   % r_RECloc(3)
end

%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
fid=fopen([result_path 'filelist.dat'],'r');
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

time0=fscanf(fid,'%8f');
status=fclose(fid);

for tcount=1:length(time0)

 zcount0=0;
 if time0(tcount)/10000000 < 1 zcount0=zcount0+1;end
 if time0(tcount)/1000000 < 1 zcount0=zcount0+1;end
 if time0(tcount)/100000 < 1 zcount0=zcount0+1;end
 if time0(tcount)/10000 < 1 zcount0=zcount0+1;end
 if time0(tcount)/1000 < 1 zcount0=zcount0+1;end
 if time0(tcount)/100 < 1 zcount0=zcount0+1;end
 if time0(tcount)/10 < 1 zcount0=zcount0+1;end
 if time0(tcount)/1 < 1 zcount0=zcount0+1;end
 time1='';
 for zcoount1=1:zcount0
  time1=[time1 '0'];
 end

 time1=[time1 int2str(time0(tcount))];
 filename=[result_path time1];
 eval(['load ' filename]);

 month0=r_time(1,2);
 month1=int2str(month0);
 date0=r_time(1,3);
 date1=int2str(date0);

%======================================================================

 obs_lat=r_h;
 obs_lat(:)=r_RECloc(1);
 obs_lon=r_h;
 obs_lon(:)=r_RECloc(2);

 degrad=pi/180.;   % conversion factor from degrees to radians
 radeg=180./pi;   % conversion factor from radians to degrees
 r_earth=6378.135; g=1.00673944;   % earth radius (km) and flatness factor

 for ii=1:length(r_h)
%  site1=[r_RECloc(1) r_RECloc(2) r_h(ii)];   % original
  site1=[site_lat site_lon r_h(ii)+site_lev];
  sinlat=sin(site1(1)*degrad);
  coslat=cos(site1(1)*degrad);
  tanlat=tan(site1(1)*degrad);
  sinlon=sin(site1(2)*degrad);
  coslon=cos(site1(2)*degrad);
  rlocgc=  [ sinlat*coslon -sinlon coslat*coslon;
             sinlat*sinlon  coslon coslat*sinlon;
            -coslat            0      sinlat    ];
  s1=r_el; s2=r_az; s3=r_range(ii);
  loc=s3*[-cos(degrad*s2)*cos(degrad*s1),...
           sin(degrad*s2)*cos(degrad*s1),...
           sin(degrad*s1)];
  gc_site1=gg2gc(site1);          % Site1 to geocentric
  gc_sp=gc_site1+(rlocgc*loc')';  % Add scattering distance in geocentric
  gg_sp=gc2gg(gc_sp);             % Transform back to geographic
  obs_lat(ii)=gg_sp(1);
  obs_lon(ii)=gg_sp(2);
 end

%======================================================================

 if tcount==1
  yy=int2str(rem(r_time(1,1),100));
  if (rem(r_time(1,1),100)) < 10 yy=['0' yy];end
  mm=int2str(r_time(1,2));
  if r_time(1,2) < 10 mm=['0' mm];end
  dd=int2str(r_time(1,3));
  if r_time(1,3) < 10 dd=['0' dd];end
  name_expr=  name_expr_set
  expr0=cap2sml(name_expr);

%======================================================================

  if (expr0 == 'cp1h') 
   if (yy =='89')|(yy == '90')|(yy =='91')
    expr0='cp1i'; name_expr='CP1I';
   end
  end
  if (expr0 == 'cp1k') 
   if (yy =='92')
    expr0='cp1j'; name_expr='CP1J';
   end
  end
  if (sdate0 == '970714') 
    expr0='enw'
  end

%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
%  output_file=[ss yy mm dd expr0 '.gdas']
  output_file=[ss yy mm dd expr0 '.gdas3']
%&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

  if exist('MCA')
   MCAs=int2str(MCA);
  else 
   MCAs='??';
  end
  if exist('MRA')
   MRAs=int2str(MRA);
  else
   MRAs='??';
  end
  MCA_com=[' Below ' MCAs ' km collision frequency is given freedom to change'];
  MRA_com=[' Above ' MRAs ' km temperture ratio is given freedom to change'];

  fid2=fopen([output_path output_file],'w');
  fprintf(fid2,'%s',name_expr);
  fprintf(fid2,'%s',' on ');
  year0=r_time(1,1);
  year1=int2str(year0);
  fprintf(fid2,'%s \n',[year1  ' ' month1 ' ' date1]);
  fprintf(fid2,'Site = ');
  fprintf(fid2,'%s \n',site0);
  fprintf(fid2,'%s \n','On fitting');
  fprintf(fid2,'%s \n',MCA_com);
  fprintf(fid2,'%s \n',MRA_com);

  if exist('veloc') & veloc == 1 
   veloc_com='Vi = positive for going away (EISCAT standard), ';
   veloc_cef=-1;
  else
   veloc_com='Vi = negative for coming closer (GUISDAP default), ';
   veloc_cef=1;
  end

  fprintf(fid2,'%s \n',[veloc_com 'lat & lon = 1e-2, Ne & col.f = log10(1e-3), O+/Ne = 1e-3']);
  title1=['date st._UT end_UT   alt  lat  lon q    Ne    Vi    Ti    Te'];
%  title1=[title1 ' Neerr Vierr Tierr Teerr    Az    El O+/Ne  colf range'];
  title1=[title1 ' Neerr Vierr Tierr Teerr    Az    El O+/Ne  colf range  Tsys TxPw'];
  title2=['MMDD HHMMSS HHMMSS    km  deg  deg     m-3   m/s     K     K'];
%  title2=[title2 '   m-3   m/s     K     K   deg   deg       rad/s    km'];
  title2=[title2 '   m-3   m/s     K     K   deg   deg       rad/s    km    K   kW '];
  fprintf(fid2,'%s \n',title1);
  fprintf(fid2,'%s \n',title2);
 end

 datef=r_time(1,2)*100+r_time(1,3);
 timef1=r_time(1,4)*10000+r_time(1,5)*100+r_time(1,6);
 timef2=r_time(2,4)*10000+r_time(2,5)*100+r_time(2,6);

 r_Te=r_param(:,2).*r_param(:,3);
 r_Teerr=r_param(:,2).*r_error(:,3)+r_param(:,3).*r_error(:,2);

 for hcount=1:length(r_h)
  altf=r_h(hcount);
  if hcount == 1, altfl=altf;end

  latf=100*obs_lat(hcount);
  lonf=100*obs_lon(hcount);
  qf=r_status(hcount);
  Nef=1e3*log10(r_param(hcount,1));
   if ~finite(Nef)|Nef>32768|Nef<-32768 Nef=-32768; end
  Vif=veloc_cef*r_param(hcount,5);
   if ~finite(Vif)|Vif>32768|Vif<-32768 Vif=-32768; end
  Tif=r_param(hcount,2);
   if ~finite(Tif)|Tif>32768|Tif<-32768 Tif=-32768; end
  Tef=r_Te(hcount);
   if ~finite(Tef)|Tef>32768|Tef<-32768 Tef=-32768; end
  Neerrf=1e3*log10(r_error(hcount,1));
   if ~finite(Neerrf)|Neerrf>32768|Neerrf<-32768 Neerrf=-32768; end
  Vierrf=r_error(hcount,5);
   if ~finite(Vierrf)|Vierrf>32768|Vierrf<-32768 Vierrf=-32768; end
  Tierrf=r_error(hcount,2);
   if ~finite(Tierrf)|Tierrf>32768|Tierrf<-32768 Tierrf=-32768; end
  Teerrf=r_Teerr(hcount);
   if ~finite(Teerrf)|Teerrf>32768|Teerrf<-32768 Teerrf=-32768; end
  ox_ratiof=1000*r_dp(hcount);
  colf=1e3*log10(r_param(hcount,4));
   if ~finite(colf) colf=-32768; end
  rangef=r_range(hcount);
  if ~finite(r_Tsys)|r_Tsys>32768|r_Tsys<-1
    r_Tsys_back=r_Tsys;
    r_Tsys=-1;
  end % error:something wrong

  fprintf(fid2,'%4.0f', datef);
  fprintf(fid2,'%7.0f', timef1);
  fprintf(fid2,'%7.0f', timef2);
  fprintf(fid2,'%6.1f', altf);
  fprintf(fid2,'%5.0f', latf);
  fprintf(fid2,'%5.0f', lonf);
  fprintf(fid2,'%2.0f', qf);
  fprintf(fid2,'%6.0f', Nef);
  fprintf(fid2,'%6.0f', Vif);
  fprintf(fid2,'%6.0f', Tif);
  fprintf(fid2,'%6.0f', Tef);
  fprintf(fid2,'%6.0f', Neerrf);
  fprintf(fid2,'%6.0f', Vierrf);
  fprintf(fid2,'%6.0f', Tierrf);
  fprintf(fid2,'%6.0f', Teerrf);
  fprintf(fid2,'%6.1f', r_az);
  fprintf(fid2,'%6.1f', r_el);
  fprintf(fid2,'%6.0f', ox_ratiof);
  fprintf(fid2,'%6.0f', colf);
%  fprintf(fid2,'%6.1f \n', rangef);
  fprintf(fid2,'%6.1f', rangef);
  fprintf(fid2,'%6.1f', r_Tsys);
  fprintf(fid2,'%5.0f \n', r_Pt/1000);
 end
 tcountmod=mod(tcount,60); % you can change 60 as you like
 if tcount == 1,
  fprintf('%6s %4s %7s %7s %7s %7s \n','tcount','date','time1','time2','alt1','alt2');end
 if tcountmod == 1,
 fprintf('%5d %5d %7d %7d %7.1f %7.1f \n',tcount,datef,timef1,timef2,altfl,altf);end

 if r_Tsys == -1,
 fprintf('%5d %5d %7d %7d %7.1f %7.1f %10.1f \n',tcount,datef,timef1,timef2,altfl,altf,r_Tsys_back)
 end
end  
 fprintf('%5d %5d %7d %7d %7.1f %7.1f \n',tcount,datef,timef1,timef2,altfl,altf);end
end 
