clear();
#Warnings3
# Komunikaty
  p1$ = "USTAW WE WŁAŚĆIWEJ POZYCJI KONSOLE\n";
  p2$ = "* TAK ABY NIE BYŁA ZASŁANIANA.\n";
  p3$ = "* NASTĘPNIE KLIKNIJ OK.";
  p$  = p1$ + p2$ + p3$;
  PopupMessage(p$);
# DEFINE PROCEDURE writeTitle: Refer to B13 & B14.
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
     printf(" analiza.sml:\n\n");
     printf("Skrypt przeznaczony do zajęć z teledetekcji.\n");
     printf("Celem jest obliczenie i analiza indeksów roślinnych\n obliczanych na podstawie wartości DN i SRFI.\n");
     printf("Jako dane wejściowe skrypt pobiera kanały:\nczerwony i podczerwony w DN i SRFI.\n");
     printf("Skrypt następnie oblicza indeks odpowiednio NDVI, IPVI, SAVI50, SAVI2\n oraz odpowiednio różnice pomiędzy wersją współczynnika obliczana z DN i SRFI\n ");
     printf("Celem ćwiczenia jest porównanie zmian w wartości\nposzególnych indeksów w zależności od rodzaju\ndanych wejściowych - DN lub SRFI\n dla powierzchni roślinnych i glebowych.\n");
     printf("Wybór powierzchni testowych jest samodzielny na podstawie kompozycji CIR lub NC (natural colors),\n.Webieramy po jednej powierzchni: roślinna - duża biomasa, glebowa - bez roślinności i mieszana - glebowo roślinna lub roślinnośc sucha i żywa.\n   ");
     printf("Należy wykonać odpowiednie opracowanie tekstowe z ilustracjami, pokazujące różnice w wartościach poszczególnych indeksów i poszczególnych rodzajów powierzchni, wskazując który indeks najabrdziej różnicuje wybrane powierzchnie ");
# deklaracja zmiennych
sleep(15)
raster R, Rc, Ir,Irc;
numeric  lin, col;
string datatype$;
#stala tla glebowego L
numeric l;
l=0.5;
# Przyporządkowanie danych zmiennym rastrowym
GetInputRaster(R);
GetInputRaster(Rc);
GetInputRaster(Ir);
GetInputRaster(Irc);
# odczytanie wymiarów i typu  rastra A w celu stworzenia 
# rastrów o tych samych wymiarach
lin = NumLins(R); col = NumCols(R); 
#stworzenie rastrów wynikowych
raster NDVIdn, NDVIsrfi;
GetOutputRaster(NDVIdn, lin, col, '32-bit float');
GetOutputRaster(NDVIsrfi, lin, col, '32-bit float');
raster IPVIdn, IPVIsrfi;
GetOutputRaster(IPVIdn, lin, col, '32-bit float');
GetOutputRaster(IPVIsrfi, lin, col, '32-bit float');
raster SAVI50dn, SAVIsrfi50;
GetOutputRaster(SAVIdn50, lin, col,'32-bit float');
GetOutputRaster(SAVIsrfi50, lin, col,'32-bit float');
raster SAVI2dn, SAVI2srfi;
GetOutputRaster(SAVI2dn, lin, col,'32-bit float');
GetOutputRaster(SAVI2srfi, lin, col,'32-bit float');
raster DIFFndvi, DIFFipvi, DIFFsavi50, DIFFsavi2;
GetOutputRaster(DIFFndvi, lin, col, "32-bit float");
GetOutputRaster(DIFFipvi, lin, col, "32-bit float");
GetOutputRaster(DIFFsavi50, lin, col, "32-bit float");
GetOutputRaster(DIFFsavi2, lin, col, "32-bit float");
#obliczenie wspolczynnikow zapisanych w formacie rastrowym

#------- ndvi --------
NDVIdn = (R-Ir)/(R+Ir);
NDVIsrfi =(Rc-Irc)/(Rc+Irc);
#------ IPVI ---------
IPVIdn = ((R-Ir)/(R+Ir)+1)/2;
IPVIsrfi =(Rc-Irc)/(Rc+Irc);
#------ SAVI50 -----
SAVIdn50=(R-Ir)*(1+l)/(R+Ir+l);
SAVIsrfi50=(Rc-Irc)*(1+l)/(Rc+Irc+l);
# ------ SAVI2--------
SAVI2dn = 0.5*(2*(Ir +1)-sqrt(pow((2*Ir+1),2)-8*(Ir-R)));
SAVI2srfi = 0.5*(2*(Irc +1)-sqrt(pow((2*Irc+1),2)-8*(Irc-Rc)));
#------- DIFF ------
DIFFndvi=NDVIdn-NDVIsrfi;
DIFFipvi=IPVIdn-IPVIsrfi;
DIFFsavi50=SAVIdn50-SAVIsrfi50;
DIFFsavi2=SAVI2dn-SAVI2srfi;
#end.