clear
start
close all


experiment=input('~~ First, give a name to your prediction run (ex: test01)=  ','s');
%experiment='test01_res2';
%disp('elisa: descomenta experiment=input')
disp(' ')
disp(['Your prediction run is:   ',experiment])

%,'coarserestest2010';


%these files contains 119000 rows and 366 column with the daily position of all virtual colonies 
%entering the Mediterranean Sea (1000 each day from 1st Jan to 30th March).

%---Define Physalia entry-time-window, years to include in database
%input("Day when Physalia started arriving at the Strait?");
disp('  ')
disp('Physalia STARTS arriving when?:')
%iyear=input('     Year for your PhysaliaSIM prediction? = ');
imonth=input('     Month of Physalia first arriving (in number)? = ');
iday=input('     Day of that month? = ');
%iyear=2010;imonth=2;iday=23;
%myyear=iyear;
%disp('elisa: descomenta input para iyear, imonth,iday')

inicio = julianday_fun(iday,imonth,2010); %Induce a poco error incluso anyo bisiesto%iyear);
clear iday,clear imonth
%inicio es 54 para 23Feb2010 iyear=input('Year of your Physalia observations? =');



%input('Day when Physalia stopped arriving at the Strait? '); %Le he quitado el 1000
disp(' ')
%disp(['In year ',int2str(iyear),', Physalia STOPS arriving on...'])
disp(['And Physalia STOPS arriving on...'])
fmonth=input(['     Month (in number)? = ']);
fday=input('     Day of that month? = ');
%fmonth=4;fday=5;
%disp('elisa: descomenta input fmonth,fday')
disp('--------')

fin = julianday_fun(fday,fmonth,2010); %Induce poco error %iyear);
clear fday,clear fmonth

%%fin=95; %5 abril2010

yearini=2008;
yearfin=2018;
%disp(' ')
%disp('elisa:  MODIFICA yearfin para incluir database hasta 2018!!!!')

disp(' ')
disp('Now, write the Lats & Lons of your Region of Interest:')
disp(' ')
minlonreg=input('Min longitude = ');
maxlonreg=input('Max longitude = ');
disp(' ')
minlatreg=input('Min latitude = ');
maxlatreg=input('Max latitude = ');
disp(' ')
%disp('elisa: descomenta los input de lat y lon, y comenta los valores abajo!!')

%Alboran
  %minlonreg=-4.7; 
  %maxlonreg=-1.9; 
  %minlatreg=36.1; 
  %maxlatreg=37.4; 

  eval(['save ',myresultsdir,'auxmat.mat experiment yearini yearfin m*reg myresultsdir'])
  eval(['save ',myresultsdir,'auxmat_',experiment,'.mat experiment yearini yearfin m*reg myresultsdir'])


disp(' ')
disp('Now you are loading the PhysaliaSIM database.')
disp('It may take some minutes ... ')
disp(' ')
disp(['Years from PhysaliaSIM database included for calculations:'])
%disp(' ')
disp(['from ',int2str(yearini),' to ',int2str(yearfin)])
disp(' ')

years=yearini:yearfin; %years you include in the database
% ---------
      
    for k=1:366 %anyo entero, aunque no entran todo el anyo!
      julday_phys((k-1)*1000+1: (k-1)*1000+1000,1)=k;
    end 
    
    %Este es el indice ind que da rango de entrada, para cada anyo
  clear ind,ind=find(julday_phys>=inicio & julday_phys<=fin);  
   julday_phys=julday_phys(ind);

%--- Only load database first time
siload=1; %input('Load/Reload database (0/1) ? = ');
   
   %----No tocar----
if siload==1
     
   %--- Build the database 
    clear data_lon*,clear data_lat*  
   for k=1:length(years)
     disp(['loading year ... ',int2str(years(k))])
     eval(['load database_modelresults_PhysaliaSIM.mat free_lon_y',int2str(years(k)),';']);
     data_longitude_free(:,:,k)=eval(['free_lon_y',int2str(years(k)),';']);
     clear free*
     eval(['load database_modelresults_PhysaliaSIM.mat free_lat_y',int2str(years(k)),';']);
     data_latitude_free(:,:,k)=eval(['free_lat_y',int2str(years(k)),';']);
     clear free*
     
   end %for k=1:length(years)
   
 
   data_longitude_free=data_longitude_free(ind,:,:);
   data_latitude_free=data_latitude_free(ind,:,:);

   disp(' ')
   disp(['Database loaded from ',int2str(years(1)),' to ',int2str(years(end))])
 
 eval(['save ',myresultsdir,'allvar_',experiment,'_database',int2str(yearini),'to',int2str(yearfin),'.mat'])
 %disp(['save allvar_',experiment,'_database',int2str(yearini),'to',int2str(yearfin),'.mat'])
else %siload=0;
 %disp(' ')
 %disp('Loading previous settings...')
 
 
 eval(['load allvar_',experiment,'_database',int2str(yearini),'to',int2str(yearfin),'.mat'])
 %disp(['load allvar_',experiment,'_database',int2str(yearini),'to',int2str(yearfin),'.mat'])
end %if siload
  
 
  %--------
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%----Compute free_physalia for many ROI
%---Define your big basin
minlonreg=-6;
maxlonreg=37;
minlatreg=29.5;
maxlatreg=46;
%--sardegna2010
%minlonreg=8.3;
%maxlonreg=10.3;
%minlatreg=38.5;
%maxlatreg=40;

clear mireg_lon,mireg_lon=[minlonreg,maxlonreg,maxlonreg,minlonreg,minlonreg];
clear mireg_lat,mireg_lat=[minlatreg,minlatreg,maxlatreg,maxlatreg,minlatreg];


disp(' ')
%disp('elisa:  MODIFICA RESOLUCION A 0.5!!!!')
intervalon=myhorizresolution; %.5;
intervalat=myhorizresolution; %.5;
disp(' ')
xmidway=round(length(minlonreg:intervalon:maxlonreg-intervalon)/2);
ymidway=round(length(minlatreg:intervalat:maxlatreg-intervalat)/2);
disp(' ')
disp('Now, the calculations using your Physalia arrival-days start ...')
i=0;
icomp=0;
for ilon=minlonreg:intervalon:maxlonreg-intervalon
  i=i+1;
  j=0;
  for jlat=minlatreg:intervalat:maxlatreg-intervalat
    icomp=icomp+1;   
     j=j+1;
    boxlon(i,j)=ilon+intervalon/2; %Punto central pixel
    boxlat(i,j)=jlat+intervalat/2;
    boxi(i,j)=i;
    boxj(i,j)=j;    
    boxid(i,j)=icomp;
    
    clear minlon,clear maxlon,clear minlat,clear maxlat
    minlon=ilon;maxlon=minlon+intervalon;
    minlat=jlat;maxlat=minlat+intervalat;
%i,j,icomp,%"pause",pause

%if i==xmidway &j==ymidway,
  %clear miroi_lon,miroi_lon=[minlon,maxlon,maxlon,minlon,minlon];
  %clear miroi_lat,miroi_lat=[minlat,minlat,maxlat,maxlat,minlat];
  %plot(mireg_lon,mireg_lat,'b','linewidth',2),hold on
  %%text((minlonreg+maxlonreg)/2,maxlatreg+maxlatreg/100,'PhysaliaSIM domain','fontsize',10,'color','b')
  %title('PhysaliaSIM domain','fontsize',10,'color','b')
  
%plot(miroi_lon,miroi_lat,'r','linewidth',2)
%text((minlon+maxlon)/2,maxlat+maxlat/100,'ROI','fontsize',10,'color','r')
%end %if i=1...

days=1:366;


% compute #arrivals each year in free_physalia
%compute the mean numeric value (for all years within the database) of arrivals to our ROI
clear free_physalia,clear mean_free,clear certainty_free_fract,clear encounter_risk
[free_physalia,mean_free,certainty_free_fract,encounter_risk]=PhysaliaSIM_mean_free_fun(minlon,maxlon,minlat,maxlat,data_longitude_free,data_latitude_free);


mat_encounter_risk(icomp,:)=encounter_risk;
mat_certainty_free_fract(icomp,:)=certainty_free_fract;


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

cualyear=1;
if i==xmidway &j==ymidway
  
    disp(' ')
    disp('Calculations are half way done ...')
    disp(' ')
       
end %if i==1...
%"pause dib",pause
%icomp
 
end %for ilon

end %for jlat


disp(' ')
disp('Calculations finished')
disp(' ')
eval(['save ',myresultsdir,'allmaps_',experiment,'_',int2str(yearini),'to',int2str(yearfin),'.mat mat* inicio fin yearini yearfin box* mireg_*'])
%disp(['save allmaps_',experiment,'_',int2str(yearini),'to',int2str(yearfin),'.mat mat* inicio fin yearini yearfin box* mireg_*'])
 

