diff --git a/Classes/AWG.m b/Classes/AWG.m
new file mode 100644
index 0000000..94c7af9
--- /dev/null
+++ b/Classes/AWG.m
@@ -0,0 +1,258 @@
+classdef AWG
+    %AWG Summary of this class goes here
+    %   Detailed explanation goes here
+
+    properties(Access=public)
+
+        kover %oversampling factor e.g. 16
+        repetitions %repeat the signal to generate a longer sequence?
+        fdac %needed
+        normalize %want to normalize at first? either 0 or 1
+        bit_resolution %bit res. of quantizer (e.g. 5 bit)
+        dac_min
+        dac_max
+
+        skew_active = 0;
+        awg_skew ; %skew vector for each output channel
+
+        lpf_active = 0;
+        lpf_type ;
+        f_cutoff ;
+        lowpass ;
+
+        signal_length;
+        H_lpf;
+
+
+    end
+
+    methods (Access=public)
+        function obj = AWG(options)
+            %AWG Construct an instance of this class
+            %   Detailed explanation goes here
+
+            arguments
+                options.preset          = [];         
+                options.kover           = 16;
+                options.repetitions     = 1;
+                options.normalize       = 1;
+                options.fdac            = 92e9;
+                options.bit_resolution  = 5.5
+                options.dac_min         = -0.5;
+                options.dac_max         = 0.5;
+                options.skew_active     = 0;
+                options.awg_skew        = 0; 
+                options.lpf_active      = 0;
+                options.lpf_type        = 0;
+                options.f_cutoff        = 32e9;
+                options.lowpass         = 1;
+            end
+
+            if isempty(options.preset)
+                obj.skew_active = 0;
+                obj.lpf_active = 0;
+
+            elseif options.preset == "M8196A"
+                % M8196A (92GBd) https://www.keysight.com/us/en/product/M8196A/92-gsa-s-arbitrary-waveform-generators.html
+                obj.dac_max         = 0.5;
+                obj.dac_min         = -.5;
+                obj.lowpass         = 1; %LP
+                obj.f_cutoff        = 32e9;
+            elseif options.preset == "M8199B"
+                %https://www.keysight.com/us/en/assets/3120-1465/data-sheets/M8199A-128-256-GSa-s-Arbitrary-Waveform-Generator.pdf
+            end
+
+            obj.kover = options.kover; %oversampling factor e.g. 16
+            obj.repetitions = options.repetitions; %repeat the signal to generate a longer sequence?
+            obj.fdac = options.fdac;
+            obj.normalize = options.normalize;%want to normalize at first? either 0 or 1
+            obj.bit_resolution = options.bit_resolution;%bit res. of quantizer (e.g. 5 bit)
+            obj.dac_min = options.dac_min;
+            obj.dac_max = options.dac_max;
+
+            obj.lpf_active = options.lpf_active;
+
+            obj.skew_active = options.skew_active;
+            obj.awg_skew = options.awg_skew;
+
+            
+
+        end
+
+        function data_out = process_channel(obj,data_in)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+            obj.signal_length = length(data_in);
+
+            if obj.normalize
+                % 5.1. Normalize the signal to 1 Vpp and set the amplitude of the signal
+                data_in = data_in/(max(data_in)-min(data_in));
+            else
+                % 5.1. Cut the Signal at -1 and 1 and scale to amplitude
+                data_in(data_in > 1) = 1;
+                data_in(data_in < -1) = -1;
+            end
+
+
+            % 1. Quantize the signal
+            if obj.bit_resolution>0
+                data_out = obj.quantization(data_in) ;
+            else
+                data_out = data_in;
+            end
+
+            % 2. Sample and hold + repeat (data_out: 1xsignal length)
+            data_out = repmat(data_out,obj.repetitions,obj.kover);
+            data_out = reshape(data_out',[],1);
+
+            % 3. Add skew
+            if obj.skew_active
+                data_out = obj.skew(data_out);
+            end         
+
+            % 4. Apply LPF on the signal
+            if obj.lpf_active
+                obj.H_lpf = obj.buildFilter(1);
+                data_out = obj.lpf(data_out) ;
+            end
+
+        end
+
+    end
+
+    methods (Access=private)
+
+        function quant_out = quantization(obj,x_in)
+
+            steps = round(2^obj.bit_resolution) + rem(round(2^obj.bit_resolution),2) ;
+
+            if isreal(x_in)
+                % shift signal and clip to quantizer intervall
+                x_in =   min(max(x_in-obj.dac_min,0),obj.dac_max-obj.dac_min);
+
+                % quantize signal
+                x_in = round((steps-1)/(obj.dac_max-obj.dac_min)*x_in);
+
+                % scale and shift back
+                quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +obj.dac_min;
+
+            else
+
+                % shift signal and clip to quantizer intervall
+                x_in = min(max(real(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min) + ...
+                    1i*min(max(imag(x_in)-obj.dac_min,0),obj.dac_max-obj.dac_min);
+
+                % quantize signal and shift back signal
+                x_in = round((steps-1)/(obj.max-obj.dac_min)*x_in);
+
+                % scale and shift back
+                quant_out = x_in*(obj.dac_max-obj.dac_min)/(steps-1) +(1+1i)*obj.dac_min;
+
+            end
+
+        end
+
+        function skew_out = skew(obj,x_in)
+
+            % Exact delay calculation
+
+            % Calculate transfer function
+            fsimu = obj.kover*obj.fdac;
+            faxis = linspace(fsimu/2, fsimu/2, length(x_in)+1);
+            faxis = fftshift(faxis(1:end-1));
+
+            skew_out = ifft(fft(x_in) .* exp(-1i*2*pi*obj.awg_skew*faxis)' );
+
+            skew_out = real(skew_out) ; % get rid of negligible imaginary part
+
+        end
+
+        function lpf_out = lpf(obj,x_in)
+            lpf_out = ifft(obj.H_lpf.*fft(x_in));
+        end
+
+        function H = buildFilter(obj,filterType)
+            
+            filtdegree = 3;
+            fsimu = obj.kover*obj.fdac;
+
+            rp = 0.5; %passband ripple
+            rs = 0.5; %stopband ripple
+
+            switch filterType
+                case 1
+
+                    % Bessel filter, impulse invariant transformed
+
+                    [B, A] = besself(filtdegree, 2*pi*obj.f_cutoff);
+                    [B ,A] = impinvar(B,A,fsimu);
+
+                case 2
+
+                    % Bessel filter, impulse bilinear transformed
+
+                    [Z, P, K] = besself(filtdegree, 2*pi*obj.f_cutoff);
+                    [Z ,P, K] = bilinear(Z,P,K,fsimu);
+                    [B ,A] = zp2tf(Z ,P ,K);
+
+                case 3
+
+                    % Butterworth filter
+                    if obj.lowpass == 1 %lowpass
+                        [B, A] = butter(filtdegree, obj.f_cutoff/(fsimu/2),'low');
+                    else % highpass
+                        [B, A] = butter(filtdegree, obj.f_cutoff/(fsimu/2),'high');
+                    end
+
+                case 4
+
+                    % Chebyshev 1 filter
+
+                    [B, A] = cheby1(filtdegree,rp, obj.f_cutoff/(fsimu/2));
+
+                case 5
+
+                    % Chebyshev 2 filter
+
+                    [B, A] = cheby2(filtdegree,rs, obj.f_cutoff/(fsimu/2));
+
+                case 6
+
+                    % Elliptic filter
+
+                    [B, A] = ellip(filtdegree,rp,rs,obj.f_cutoff/(fsimu/2));
+
+                case 7
+
+                    % Hamming filter
+
+                    g=(filtdegree-1)/2;
+                    wc=obj.f_cutoff/(fsimu/2);
+                    B = wc*sinc(wc*(-g:g)).*hamming(filtdegree)';
+                    A=1;
+
+                case 8
+
+                    % Raised Cosine filter
+
+                    B = firrcos(filtdegree,obj.f_cutoff,para.df,fsimu);
+                    A=1;
+
+                case 9
+
+                    % Sinc filter
+
+                    g=(filtdegree-1)/2;
+                    wc=obj.f_cutoff/(fsimu/2);
+                    B = wc*sinc(wc*(-g:g));
+                    A=1;
+
+            end
+
+             H = freqz(B, A, obj.repetitions*obj.kover*obj.signal_length,'whole');
+
+
+        end
+
+    end
+end
diff --git a/Classes/Amplifier.m b/Classes/Amplifier.m
new file mode 100644
index 0000000..d11eaee
--- /dev/null
+++ b/Classes/Amplifier.m
@@ -0,0 +1,140 @@
+classdef Amplifier
+    %UNTITLED Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties
+        type 
+        amp_mode
+        amplification_db
+
+        nase_mode
+        noifig
+
+        saturation_mode
+        saturation_power
+
+        fsimu
+
+    end
+    
+    methods
+        function obj = Amplifier(options)
+            %UNTITLED Construct an instance of this class
+            %   Detailed explanation goes here
+            arguments
+                options.type
+                options.amp_mode
+                options.amplification_db
+
+                options.nase_mode = 0;
+                options.noifig = 0;
+
+                options.saturation_mode = 0;
+                options.saturation_power = 0;
+
+                options.fsimu = []
+            end
+
+            obj.type = options.type;
+            obj.amp_mode = options.amp_mode;
+            obj.amplification_db = options.amplification_db;
+
+            obj.nase_mode = options.nase_mode;
+            obj.noifig = options.noifig;
+
+            obj.saturation_mode = options.saturation_mode;
+            obj.saturation_power = options.saturation_power;
+
+            obj.fsimu = options.fsimu;
+
+
+            obj.saturation_power=1/1000*10^(obj.saturation_power/10);
+
+
+        end
+        
+        function [y_out,nase] = process(obj,x_in,optional)
+
+            arguments
+                obj
+                x_in = [];
+
+                optional.nase = 0;
+            end
+
+            %calc gain for output power mode, amp mode and saturated mode
+            a_lin = obj.calculateGain(x_in);
+
+            y_out = a_lin * x_in;
+
+            if obj.type == "ideal"
+
+                %don't add/ remove noise, but scale it accordingly
+                if optional.nase ~= 0
+                    nase=state.a^2*optional.nase;
+                end
+
+            elseif obj.type == "edfa"
+
+                %calculate ASE-noise
+                if optional.nase ~= 0
+                    Nase_old = state.gain^2*optional.nase ;
+
+                    h = Constant.Planck;
+                    c = Constant.LightSpeed;
+
+                    Nase_new = 0.5* 10^(obj.noifig/10) * h*c/(lambda_T*1e-9) * (a_lin^2-1) ;
+
+                    nase = Nase_old + Nase_new ; 
+                end
+
+            end
+
+            pow_in_lin =  mean(abs(x_in.^2)) ; 
+            pow_in_dbm = 10*log10(pow_in_lin)+30;
+            pow_out_lin = mean(abs(y_out.^2)) ; 
+            pow_out_dbm = 10*log10(pow_out_lin)+30;
+
+            if obj.amp_mode == "output_power"
+                seemsright = pow_out_dbm == obj.amplification_db;
+            else obj.amp_mode == "gain"
+                seemsright = pow_out_dbm == pow_in_dbm+ obj.amplification_db;
+            end
+
+            if ~seemsright 
+                warning("Amplifier output not correct, please check the reason");
+            end
+
+        end
+
+        function a_lin = calculateGain(obj,xin)
+
+            if obj.amp_mode == "output_power"
+
+                %get linear gain for output power mode
+                pow_in = mean(abs(xin.^2)) ; % lin input power
+                pow_out = 10^(obj.amplification_db/10 - 3) ; % dBm to lin
+                a_lin = sqrt(pow_out/pow_in) ;
+
+            elseif obj.amp_mode == "gain"
+
+                %get linear gain for classic gain mode 
+                a_lin=10^(obj.amplification_db/20);
+                
+            end
+
+            if obj.saturation_mode 
+                Pin = sum(mean((abs(xin)).^2, 2)) ;
+
+                gain_power=fzero(inline(['log(G/' num2str(a_lin^2,'%1.16e') ')/(1-G)-log(2)*' num2str(Pin,'%1.16e') '/' ...
+                num2str(obj.saturation_power,'%1.16e')],'G'),[1+eps,a_lin^2]);
+
+            end
+
+        end
+
+
+
+    end
+end
+
diff --git a/Classes/EML.m b/Classes/EML.m
new file mode 100644
index 0000000..cfe4efd
--- /dev/null
+++ b/Classes/EML.m
@@ -0,0 +1,135 @@
+classdef EML
+    %EML External Modulated Laser
+    %   Detailed explanation goes here
+    
+    properties
+        %class input objms
+        mode
+        fsimu %in Hz
+        lambda %in nm
+        power %in dBm
+        linewidth %in Hz
+        ampl_imbal 
+        pha_imbal
+        bias 
+        u_pi
+
+        %on instance creation
+        field
+        noisefactor
+        phase %eigentlich nur wichtig wenn man mehrere blöcke durchsimuliert...
+        real_factor
+        imag_factor
+
+        %during simulation
+        signal_len
+
+    end
+    
+    methods
+        function obj = EML(options)
+            %EML Construct an instance of this class
+            %   Detailed explanation goes here
+            arguments
+                options.mode;
+                options.fsimu;
+                options.lambda;
+                options.power;
+                options.linewidth = 0;
+                options.ampl_imbal = 0;
+                options.pha_imbal = 0;
+                options.bias;
+                options.u_pi;
+            end
+            
+            obj.mode = options.mode;
+            obj.fsimu = options.fsimu;
+            obj.lambda = options.lambda;
+            obj.power = options.power;
+            obj.linewidth = options.linewidth;
+            obj.ampl_imbal = options.ampl_imbal;
+            obj.pha_imbal = options.pha_imbal;
+            obj.bias = options.bias;
+            obj.u_pi = options.u_pi;
+
+            obj.field=sqrt(10^(obj.power/10-3)); %dbm to sqrt(mw)
+            
+            obj.noisefactor=sqrt(2*pi*obj.linewidth/obj.fsimu);
+            obj.phase=0;
+
+            if obj.mode==4 || obj.mode==5
+                obj.real_factor=(1+obj.ampl_imbal/2)*exp(1i*obj.pha_imbal/2);
+                obj.imag_factor=(1-obj.ampl_imbal/2)*exp(1i*(pi/2-obj.pha_imbal/2));
+            end
+
+        end
+        
+        function [opt_out,obj] = process(obj,elec_in)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+
+            %needed later
+            obj.signal_len = length(elec_in);
+
+            % generate phase noise of laser
+            % scale with signal energy in sqrt(mw)
+            if obj.linewidth ~= 0
+                ph_noi = obj.createPhaseNoise;
+                laserfield = obj.field.*exp(1i*ph_noi);
+            else
+                laserfield = obj.field;
+            end
+            
+            %modulate the laserfield with the electrical signal
+            opt_out = obj.externalmodulation(laserfield,elec_in);
+
+            %remember phase (! you need to receive the altered eml object in you sim program !)
+            obj.phase = ph_noi(end); 
+
+        end
+
+        function noi = createPhaseNoise(obj)
+            %create random vector
+            noi = randn(obj.signal_len,1);
+            %scale with noisefactor
+            noi = noi * obj.noisefactor;
+            %cumsum to accumulate noise over time vector
+            noi = cumsum(noi);
+            %add phase from previous block/ loop of simulation
+            noi = noi + obj.phase;
+        end
+
+        function modulated_laserfield = externalmodulation(obj,laserfield,electrical_in)
+
+            switch obj.mode
+
+                case   1 % linear IM
+                    modulated_laserfield = laserfield.*sqrt(real(electrical_in + obj.bias)/obj.u_pi);
+
+                case 2   % IM with hMZM function
+                    modulated_laserfield = laserfield.*cos(pi/2*(real(electrical_in)+obj.bias)/obj.u_pi);
+
+                case 3 % PMsi
+                    modulated_laserfield = laserfield.*exp(1i*pi*real(electrical_in)/obj.u_pi);
+
+                case 4 % linear IQ-Modulator
+                    if obj.ampl_imbal ==0 && obj.pha_imbal  == 0
+                        % No IQ Imbalance
+                        modulated_laserfield=laserfield.*electrical_in/obj.u_pi;
+                    else
+                        % with IQ Imbalance
+                        modulated_laserfield = laserfield.*(real(electrical_in)*obj.real_factor/obj.u_pi...
+                            +imag(electrical_in)*obj.imag_factor/obj.u_pi);
+                    end
+
+                case 5 % IQ-Modulator with cosine function
+                    modulated_laserfield = laserfield/sqrt(2).*(cos(pi/2*(real(electrical_in)+obj.bias)/obj.u_pi)*obj.real_factor...
+                        +cos(pi/2*(imag(electrical_in)+obj.bias)/obj.u_pi)*obj.imag_factor);
+
+            end
+
+        end
+
+    end
+end
+
diff --git a/Classes/Electricalsignal.m b/Classes/Electricalsignal.m
new file mode 100644
index 0000000..b6f2b77
--- /dev/null
+++ b/Classes/Electricalsignal.m
@@ -0,0 +1,25 @@
+classdef Electricalsignal < Signal
+    %ELECTRICALSIGNAL Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties
+        
+    end
+    
+    methods
+        function obj = Electricalsignal(signal, fsym, fsimu)
+            %ELECTRICALSIGNAL Construct an instance of this class
+            %   Detailed explanation goes here
+            obj = obj@Signal(signal, fsym, fsimu);
+            
+        end
+        
+        function outputArg = method1(obj,inputArg)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+            outputArg = obj.Property1 + inputArg;
+        end
+
+    end
+end
+
diff --git a/Classes/Fiber.m b/Classes/Fiber.m
new file mode 100644
index 0000000..aca6995
--- /dev/null
+++ b/Classes/Fiber.m
@@ -0,0 +1,51 @@
+classdef Fiber
+    %FIBER Summary of this class goes here
+    %   Detailed explanation goes here
+
+    properties
+        fiber_length
+        alpha
+        D
+        Dslope
+        lambda0
+        gamma
+        dphimax
+    end
+
+    methods
+        function obj = Fiber(options)
+            %FIBER Construct an instance of this class
+            %   Detailed explanation goes here
+            arguments
+                options.fiber_length = 0
+                options.alpha = 0.2
+                options.D = 17
+                options.Dslope = 0.06
+                options.lambda0 = 1550
+                options.gamma = 0.0013
+                options.dphimax = 5e-3
+            end
+
+            obj.fiber_length = options.fiber_length*1000; %km
+            obj.alpha = options.alpha;
+            obj.D =options.D*1e-6;
+            obj.Dslope =options.Dslope;
+            obj.lambda0 =options.lambda0;
+            obj.gamma =options.gamma;
+            obj.dphimax =options.dphimax;
+
+        end
+
+        function opt_out = process(obj,opt_in)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+
+            %attenuate nase
+
+
+
+
+        end
+    end
+end
+
diff --git a/Classes/Filter.m b/Classes/Filter.m
new file mode 100644
index 0000000..f5e82af
--- /dev/null
+++ b/Classes/Filter.m
@@ -0,0 +1,139 @@
+classdef Filter
+    %FILTER Summary of this class goes here
+    %   Detailed explanation goes here
+
+    properties
+        H
+        filterType
+        f_cutoff
+
+
+        signal_length
+        filtdegree
+        passband_ripple 
+        stopband_ripple
+        fdac
+    end
+
+    methods
+        function obj = Filter(options)
+            %FILTER Construct an instance of this class
+            %   Detailed explanation goes here
+            arguments
+                options.filterType = 1;
+                options.f_cutoff   = 0;
+                options.fdac       = 0;
+                options.filtdegree = 3;
+                options.passband_ripple = 0.5;
+                options.stopband_ripple = 0.5;
+            end
+            
+
+            obj.filterType = options.filterType;
+
+
+            obj.f_cutoff        = options.f_cutoff;
+            obj.filtdegree      = options.filtdegree;
+            obj.passband_ripple = options.passband_ripple;
+            obj.stopband_ripple = options.stopband_ripple;
+            obj.fdac            = options.fdac;
+
+        end
+
+        function yout = process(obj,xin)
+
+            obj.signal_length = length(xin);
+
+            obj.H = obj.buildFilter(obj.filterType);
+
+            yout = obj.applyFilter(xin);
+
+        end
+
+        function y_filtered = applyFilter(obj,xin)
+
+            y_filtered = ifft(obj.H.*fft(xin));
+
+        end
+
+        function H = buildFilter(obj,filterType)
+
+            rp = obj.passband_ripple; %passband ripple
+            rs = obj.stopband_ripple; %stopband ripple
+
+            switch filterType
+                case 1
+
+                    % Bessel filter, impulse invariant transformed
+
+                    [B, A] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
+                    [B ,A] = impinvar(B,A,obj.fdac);
+
+                case 2
+
+                    % Bessel filter, impulse bilinear transformed
+
+                    [Z, P, K] = besself(obj.filtdegree, 2*pi*obj.f_cutoff);
+                    [Z ,P, K] = bilinear(Z,P,K,obj.fdac);
+                    [B ,A] = zp2tf(Z ,P ,K);
+
+                case 3
+
+                    % Butterworth filter
+                    if obj.lowpass == 1 %lowpass
+                        [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'low');
+                    else % highpass
+                        [B, A] = butter(obj.filtdegree, obj.f_cutoff/(obj.fdac/2),'high');
+                    end
+
+                case 4
+
+                    % Chebyshev 1 filter
+
+                    [B, A] = cheby1(obj.filtdegree,rp, obj.f_cutoff/(obj.fdac/2));
+
+                case 5
+
+                    % Chebyshev 2 filter
+
+                    [B, A] = cheby2(obj.filtdegree,rs, obj.f_cutoff/(obj.fdac/2));
+
+                case 6
+
+                    % Elliptic filter
+
+                    [B, A] = ellip(obj.filtdegree,rp,rs,obj.f_cutoff/(obj.fdac/2));
+
+                case 7
+
+                    % Hamming filter
+
+                    g=(obj.filtdegree-1)/2;
+                    wc=obj.f_cutoff/(obj.fdac/2);
+                    B = wc*sinc(wc*(-g:g)).*hamming(obj.filtdegree)';
+                    A=1;
+
+                case 8
+
+                    % Raised Cosine filter
+
+                    B = firrcos(obj.filtdegree,obj.f_cutoff,para.df,obj.fdac);
+                    A=1;
+
+                case 9
+
+                    % Sinc filter
+
+                    g=(obj.filtdegree-1)/2;
+                    wc=obj.f_cutoff/(obj.fdac/2);
+                    B = wc*sinc(wc*(-g:g));
+                    A=1;
+
+            end
+
+            H = freqz(B, A, obj.signal_length,'whole');
+
+        end
+    end
+end
+
diff --git a/Classes/Opticalsignal.m b/Classes/Opticalsignal.m
new file mode 100644
index 0000000..cf2f648
--- /dev/null
+++ b/Classes/Opticalsignal.m
@@ -0,0 +1,28 @@
+classdef Opticalsignal < Signal
+    %OPTICALSIGNAL Summary of this class goes here
+    %   Detailed explanation goes here
+
+    properties
+        nase
+        lambda_nm
+    end
+
+    methods
+        function obj = Opticalsignal(signal, fsym, fsimu)
+            %OPTICALSIGNAL Construct an instance of this class
+            %   Detailed explanation goes here
+            obj = obj@Signal(signal, fsym, fsimu);
+
+            obj.nase = 0;
+            obj.lambda_nm = [];
+            
+        end
+
+        function outputArg = method1(obj,inputArg)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+            outputArg = obj.Property1 + inputArg;
+        end
+    end
+end
+
diff --git a/Classes/Signal.m b/Classes/Signal.m
new file mode 100644
index 0000000..3c6c202
--- /dev/null
+++ b/Classes/Signal.m
@@ -0,0 +1,48 @@
+classdef Signal
+    %SIGNAL Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties
+        signal
+        fsym
+        fsimu
+    end
+    
+    methods
+        function obj = Signal(signal, fsym, fsimu)
+            %SIGNAL Construct an instance of this class
+            %   Detailed explanation goes here
+            obj.signal = signal;
+            
+            obj.fsym = fsym;
+            obj.fsimu = fsimu;
+
+        end
+
+        function [e_sig, nase, lambda_nm] = Electricalsignal(obj)
+
+            if isa(obj,'Opticalsignal')
+                %convert to electrical
+                disp("Convert signal: opt. -> elec. Fetch all properties (e.g. nase)!");
+                nase = obj.nase;
+                lambda_nm = obj.lambda_nm;
+                e_sig = Electricalsignal(obj.signal,obj.fsym, obj.fsimu);
+            end
+            
+        end
+
+        function o_sig = Opticalsignal(obj)
+            if isa(obj,'Electricalsignal')
+                %convert to optical
+                disp("Convert signal: elec. -> opt.!");
+                o_sig = Opticalsignal(obj.signal,obj.fsym, obj.fsimu);
+            end
+        end
+
+        function outputArg = method1(obj,inputArg)
+
+        end
+
+    end
+end
+
diff --git a/Datatypes/emlmodes.m b/Datatypes/emlmodes.m
new file mode 100644
index 0000000..9486815
--- /dev/null
+++ b/Datatypes/emlmodes.m
@@ -0,0 +1,11 @@
+classdef emlmodes < int32
+
+   enumeration
+      im_linear  (1)
+      im_cosinus (2)
+      pm         (3)
+      iq_linear  (4)
+      iq_cosinus (5)
+   end
+   
+end
\ No newline at end of file
diff --git a/Datatypes/filtertypes.m b/Datatypes/filtertypes.m
new file mode 100644
index 0000000..3633519
--- /dev/null
+++ b/Datatypes/filtertypes.m
@@ -0,0 +1,15 @@
+classdef filtertypes < int32
+
+   enumeration
+      bessel_inp  (1)
+      bessel_bilin(2)
+      butterworth (3)
+      chebyshev1  (4)
+      chebyshev2  (5)
+      elliptic    (6)
+      hamming     (7)
+      racos       (8)
+      sinc        (9)
+   end
+
+end
diff --git a/Fiber.asv b/Fiber.asv
new file mode 100644
index 0000000..23c8d17
--- /dev/null
+++ b/Fiber.asv
@@ -0,0 +1,52 @@
+classdef Fiber
+    %FIBER Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties
+            fiber_length
+            alpha
+            D
+            Dslope
+            lambda0
+            gamma
+            dphimax
+    end
+    
+    methods
+        function obj = Fiber(options)
+            %FIBER Construct an instance of this class
+            %   Detailed explanation goes here
+            arguments
+                options.fiber_length = 0
+                options.alpha = 0.2
+                options.D = 17
+                options.Dslope = 0.06
+                options.lambda0 = 1550
+                options.gamma = 0.0013
+                options.dphimax = 5e-3
+            end
+
+            obj.fiber_length
+            obj.alpha
+            obj.D
+            Dslope
+            lambda0
+            gamma
+            dphimax
+
+
+        end
+        
+        function opt_out = process(obj,opt_in)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+            
+            %attenuate nase
+
+
+
+
+        end
+    end
+end
+
diff --git a/Signal.asv b/Signal.asv
new file mode 100644
index 0000000..5fa6735
--- /dev/null
+++ b/Signal.asv
@@ -0,0 +1,31 @@
+classdef Signal
+    %SIGNAL Summary of this class goes here
+    %   Detailed explanation goes here
+    
+    properties
+        signal
+        nase
+        fsym
+        fsimu
+        lambda
+    end
+    
+    methods
+        function obj = Signal(signal, fsym, fsimu)
+            %SIGNAL Construct an instance of this class
+            %   Detailed explanation goes here
+            obj.signal = signal;
+            obj.nase = 0;
+            obj.fsym = fsym;
+            obj.fsimu = fsimu;
+            obj.lambda = 
+        end
+        
+        function outputArg = method1(obj,inputArg)
+            %METHOD1 Summary of this method goes here
+            %   Detailed explanation goes here
+            outputArg = obj.Property1 + inputArg;
+        end
+    end
+end
+
diff --git a/comm_tb.asv b/comm_tb.asv
new file mode 100644
index 0000000..e7c829d
--- /dev/null
+++ b/comm_tb.asv
@@ -0,0 +1,146 @@
+O = 10; %order of prbs
+N = 64; %length of prbs
+[~,seed] = prbs(O,1); %initialize first seed of prbs
+
+% Modulation
+M = 4;
+bitpattern = zeros(N,log2(M));
+
+% Symbol Rate
+fsym = 56e9;
+% DAC Rate
+fdac = 120e9;
+% Simulation oversampling rate "k";
+kover = 16;
+% Simulation frequency in "analog domain"
+fsimu = kover * fdac ;
+
+
+
+%SIMULATE
+
+for i = 1:log2(M)
+    [bitpattern(:,i),seed] = prbs(O,N,seed);
+end
+
+pamData     = pam_mapping(bitpattern,M,0);
+
+shapedData  = applyPulseShaping(pamData,fsym,fdac);
+
+
+
+awg = AWG('preset','M8196A','fdac',fdac,'kover',kover,'lpf_active',1);
+awgSignal   = awg.process_channel(shapedData);
+
+fil = Filter('filtdegree',4,"f_cutoff",50e9,"fdac",fdac,"filterType",filtertypes.bessel_bilin);
+filtered = fil.process(awgSignal);
+
+u_pi = 3.5;
+vbias = (0.5*u_pi)-u_pi;
+eml = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000);
+laserfield = eml.process(filtered); 
+
+att = Amplifier("amplification_db",10,"amp_mode","gain","type","ideal","saturation_mode",0,'saturation_power',10);
+att_out = att.process(laserfield);
+
+%fib = Fiber("fiber_length",2,"alpha",0.2,"D",17,"lambda0",1550);
+
+
+figure;
+plot(shapedData);
+hold on
+plot(awgSignal,'DisplayName','skew 0');
+plot(filtered,'DisplayName','filtered');
+plot(abs(laserfield),'DisplayName','laser');
+plot(abs(att_out),'DisplayName','att_out');
+
+hold off
+
+
+function pam_sig = pam_mapping(bitpattern, M, unipolar)
+
+switch log2(M)
+    case 1
+        % 2-ASK: BPSK / OOK
+        pam_sig=bitpattern(:,1);
+
+        if unipolar==0
+            pam_sig=2*pam_sig-1;
+        end
+
+    case 2
+        % 4-ASK:
+        pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2));
+
+        if unipolar==0
+            pam_sig=2*pam_sig-3;
+        end
+
+
+    case 3
+        % 8-ASK:
+        x1 = bitpattern(:,1);
+        x2 = (bitpattern(:,1)==bitpattern(:,3));
+        x3 = x2~=bitpattern(:,2);
+
+        pam_sig = 4*x1 + 2*x2 + x3;
+
+        if unipolar==0
+            pam_sig=2*pam_sig-7;
+        end
+
+
+    case 4
+        % 16-ASK:
+        x1 = bitpattern(:,1);
+        x2 = (bitpattern(:,1)==bitpattern(:,4));
+        x3 = x2~=bitpattern(:,3);
+        x4 = x3~=bitpattern(:,2);
+
+        pam_sig = 8*x1 + 4*x2 + 2*x3 + x4;
+
+        if unipolar==0
+            pam_sig=2*pam_sig-15;
+        end
+
+end
+end
+
+function yout = applyPulseShaping(xin,fsym,fdac)
+
+    if ~rem(fdac,fsym)
+        %ist ein Vielfaches
+        sps = fdac / fsym;
+        up = sps;
+        dn = 1;
+    else
+        %ist kein Vielfaches
+        up = fdac / gcd(fdac, fsym);
+        dn = fsym / gcd(fdac, fsym);
+        sps= up;
+    end
+
+    %Bau das Filter (hier rrc)
+    racos_len = 2048;
+    alpha = 0.1;
+    h = rcosdesign(alpha,racos_len,sps);
+    
+    %Apply Filter using Matlab build in fctn. 
+    yout = upfirdn(xin,h,up,dn);
+    
+    %cut signal, which is longer due to fir filter 
+    st = up/dn*racos_len/2; %we need to cut y_out 
+    en = st + (length(xin)*up/dn) -1;
+
+    yout = yout(st:en);
+
+    %Check output integrity 
+    if (up/dn * length(xin)) ~= length(yout)
+        warning('Check signal length after pulse shaping');
+    end
+
+end
+
+function yout = applyBandwidthLimitation(xin)
+    data_out=ifft(repmat(state.H,1,size(data_in,1)).*fft(data_in.')).';
+end
diff --git a/comm_tb.m b/comm_tb.m
new file mode 100644
index 0000000..785399e
--- /dev/null
+++ b/comm_tb.m
@@ -0,0 +1,146 @@
+O = 10; %order of prbs
+N = 64; %length of prbs
+[~,seed] = prbs(O,1); %initialize first seed of prbs
+
+% Modulation
+M = 4;
+bitpattern = zeros(N,log2(M));
+
+% Symbol Rate
+fsym = 56e9;
+% DAC Rate
+fdac = 120e9;
+% Simulation oversampling rate "k";
+kover = 16;
+% Simulation frequency in "analog domain"
+fsimu = kover * fdac ;
+
+
+
+%SIMULATE
+
+for i = 1:log2(M)
+    [bitpattern(:,i),seed] = prbs(O,N,seed);
+end
+
+pamData     = pam_mapping(bitpattern,M,0);
+
+shapedData  = applyPulseShaping(pamData,fsym,fdac);
+
+
+
+awg = AWG('preset','M8196A','fdac',fdac,'kover',kover,'lpf_active',1);
+awgSignal   = awg.process_channel(shapedData);
+
+fil = Filter('filtdegree',4,"f_cutoff",50e9,"fdac",fdac,"filterType",filtertypes.bessel_bilin);
+filtered = fil.process(awgSignal);
+
+u_pi = 3.5;
+vbias = (0.5*u_pi)-u_pi;
+eml = EML("mode",emlmodes.im_cosinus,"power",10,"fsimu",fsimu,"lambda",1550,"bias",vbias,"u_pi",u_pi,"linewidth",1000000);
+laserfield = eml.process(filtered); 
+
+att = Amplifier("amplification_db",10,"amp_mode","gain","type","ideal","saturation_mode",0,'saturation_power',10);
+att_out = att.process(laserfield);
+
+%fib = Fiber("fiber_length",2,"alpha",0.2,"D",17,"lambda0",1550);
+
+
+figure;
+plot(shapedData);
+hold on
+plot(awgSignal,'DisplayName','skew 0');
+plot(filtered,'DisplayName','filtered');
+plot(abs(laserfield),'DisplayName','laser');
+plot(abs(att_out),'DisplayName','att_out');
+
+hold off
+
+
+function pam_sig = pam_mapping(bitpattern, M, unipolar)
+
+switch log2(M)
+    case 1
+        % 2-ASK: BPSK / OOK
+        pam_sig=bitpattern(:,1);
+
+        if unipolar==0
+            pam_sig=2*pam_sig-1;
+        end
+
+    case 2
+        % 4-ASK:
+        pam_sig=2*bitpattern(:,1)+(bitpattern(:,1)==bitpattern(:,2));
+
+        if unipolar==0
+            pam_sig=2*pam_sig-3;
+        end
+
+
+    case 3
+        % 8-ASK:
+        x1 = bitpattern(:,1);
+        x2 = (bitpattern(:,1)==bitpattern(:,3));
+        x3 = x2~=bitpattern(:,2);
+
+        pam_sig = 4*x1 + 2*x2 + x3;
+
+        if unipolar==0
+            pam_sig=2*pam_sig-7;
+        end
+
+
+    case 4
+        % 16-ASK:
+        x1 = bitpattern(:,1);
+        x2 = (bitpattern(:,1)==bitpattern(:,4));
+        x3 = x2~=bitpattern(:,3);
+        x4 = x3~=bitpattern(:,2);
+
+        pam_sig = 8*x1 + 4*x2 + 2*x3 + x4;
+
+        if unipolar==0
+            pam_sig=2*pam_sig-15;
+        end
+
+end
+end
+
+function yout = applyPulseShaping(xin,fsym,fdac)
+
+    if ~rem(fdac,fsym)
+        %ist ein Vielfaches
+        sps = fdac / fsym;
+        up = sps;
+        dn = 1;
+    else
+        %ist kein Vielfaches
+        up = fdac / gcd(fdac, fsym);
+        dn = fsym / gcd(fdac, fsym);
+        sps= up;
+    end
+
+    %Bau das Filter (hier rrc)
+    racos_len = 2048;
+    alpha = 0.1;
+    h = rcosdesign(alpha,racos_len,sps);
+    
+    %Apply Filter using Matlab build in fctn. 
+    yout = upfirdn(xin,h,up,dn);
+    
+    %cut signal, which is longer due to fir filter 
+    st = up/dn*racos_len/2; %we need to cut y_out 
+    en = st + (length(xin)*up/dn) -1;
+
+    yout = yout(st:en);
+
+    %Check output integrity 
+    if (up/dn * length(xin)) ~= length(yout)
+        warning('Check signal length after pulse shaping');
+    end
+
+end
+
+function yout = applyBandwidthLimitation(xin)
+    data_out=ifft(repmat(state.H,1,size(data_in,1)).*fft(data_in.')).';
+end
diff --git a/physconst/+Constant/+HA/Boltzmann.m b/physconst/+Constant/+HA/Boltzmann.m
new file mode 100644
index 0000000..88fa71a
--- /dev/null
+++ b/physconst/+Constant/+HA/Boltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=Boltzmann()
+
+v=1;
+u=9.1e-7;
+
+end
diff --git a/physconst/+Constant/+HA/ConductanceQuantum.m b/physconst/+Constant/+HA/ConductanceQuantum.m
new file mode 100644
index 0000000..6b566ab
--- /dev/null
+++ b/physconst/+Constant/+HA/ConductanceQuantum.m
@@ -0,0 +1,6 @@
+function [v,u]=ConductanceQuantum()
+
+v=1/pi;
+u=1e-10;
+
+end
diff --git a/physconst/+Constant/+HA/Coulomb.m b/physconst/+Constant/+HA/Coulomb.m
new file mode 100644
index 0000000..983cead
--- /dev/null
+++ b/physconst/+Constant/+HA/Coulomb.m
@@ -0,0 +1,6 @@
+function [v,u]=Coulomb()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+HA/ElectronMass.m b/physconst/+Constant/+HA/ElectronMass.m
new file mode 100644
index 0000000..401a531
--- /dev/null
+++ b/physconst/+Constant/+HA/ElectronMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ElectronMass()
+
+v=1;
+u=4.4e-8;
+
+end
diff --git a/physconst/+Constant/+HA/ElementaryCharge.m b/physconst/+Constant/+HA/ElementaryCharge.m
new file mode 100644
index 0000000..432a486
--- /dev/null
+++ b/physconst/+Constant/+HA/ElementaryCharge.m
@@ -0,0 +1,6 @@
+function [v,u]=ElementaryCharge()
+
+v=1;
+u=2.2e-8;
+
+end
diff --git a/physconst/+Constant/+HA/FineStructure.m b/physconst/+Constant/+HA/FineStructure.m
new file mode 100644
index 0000000..a75b283
--- /dev/null
+++ b/physconst/+Constant/+HA/FineStructure.m
@@ -0,0 +1,6 @@
+function [v,u]=FineStructure()
+
+v=7.2973525698e-3;
+u=2.4e-12;
+
+end
diff --git a/physconst/+Constant/+HA/FluxQuantum.m b/physconst/+Constant/+HA/FluxQuantum.m
new file mode 100644
index 0000000..ba3ae89
--- /dev/null
+++ b/physconst/+Constant/+HA/FluxQuantum.m
@@ -0,0 +1,6 @@
+function [v,u]=FluxQuantum()
+
+v=pi;
+u=6.9e-8;
+
+end
diff --git a/physconst/+Constant/+HA/LightSpeed.m b/physconst/+Constant/+HA/LightSpeed.m
new file mode 100644
index 0000000..24d645b
--- /dev/null
+++ b/physconst/+Constant/+HA/LightSpeed.m
@@ -0,0 +1,6 @@
+function [v,u]=LightSpeed()
+
+v=137.03599907;
+u=4.4e-8;
+
+end
diff --git a/physconst/+Constant/+HA/Planck.m b/physconst/+Constant/+HA/Planck.m
new file mode 100644
index 0000000..685f195
--- /dev/null
+++ b/physconst/+Constant/+HA/Planck.m
@@ -0,0 +1,6 @@
+function [v,u]=Planck()
+
+v=2*pi;
+u=7.5e-8;
+
+end
diff --git a/physconst/+Constant/+HA/ProtonMass.m b/physconst/+Constant/+HA/ProtonMass.m
new file mode 100644
index 0000000..7a336e2
--- /dev/null
+++ b/physconst/+Constant/+HA/ProtonMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ProtonMass()
+
+v=1836.15267;
+u=1.6e-4;
+
+end
diff --git a/physconst/+Constant/+HA/ReducedPlanck.m b/physconst/+Constant/+HA/ReducedPlanck.m
new file mode 100644
index 0000000..6792f15
--- /dev/null
+++ b/physconst/+Constant/+HA/ReducedPlanck.m
@@ -0,0 +1,6 @@
+function [v,u]=ReducedPlanck()
+
+v=1;
+u=1.2e-8;
+
+end
diff --git a/physconst/+Constant/+P/Boltzmann.m b/physconst/+Constant/+P/Boltzmann.m
new file mode 100644
index 0000000..88fa71a
--- /dev/null
+++ b/physconst/+Constant/+P/Boltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=Boltzmann()
+
+v=1;
+u=9.1e-7;
+
+end
diff --git a/physconst/+Constant/+P/Coulomb.m b/physconst/+Constant/+P/Coulomb.m
new file mode 100644
index 0000000..983cead
--- /dev/null
+++ b/physconst/+Constant/+P/Coulomb.m
@@ -0,0 +1,6 @@
+function [v,u]=Coulomb()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+P/Gravitational.m b/physconst/+Constant/+P/Gravitational.m
new file mode 100644
index 0000000..2e47506
--- /dev/null
+++ b/physconst/+Constant/+P/Gravitational.m
@@ -0,0 +1,6 @@
+function [v,u]=Gravitational()
+
+v=1;
+u=4.7e-5;
+
+end
diff --git a/physconst/+Constant/+P/LightSpeed.m b/physconst/+Constant/+P/LightSpeed.m
new file mode 100644
index 0000000..02df4a3
--- /dev/null
+++ b/physconst/+Constant/+P/LightSpeed.m
@@ -0,0 +1,6 @@
+function [v,u]=LightSpeed()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+P/Planck.m b/physconst/+Constant/+P/Planck.m
new file mode 100644
index 0000000..685f195
--- /dev/null
+++ b/physconst/+Constant/+P/Planck.m
@@ -0,0 +1,6 @@
+function [v,u]=Planck()
+
+v=2*pi;
+u=7.5e-8;
+
+end
diff --git a/physconst/+Constant/+P/ReducedPlanck.m b/physconst/+Constant/+P/ReducedPlanck.m
new file mode 100644
index 0000000..6792f15
--- /dev/null
+++ b/physconst/+Constant/+P/ReducedPlanck.m
@@ -0,0 +1,6 @@
+function [v,u]=ReducedPlanck()
+
+v=1;
+u=1.2e-8;
+
+end
diff --git a/physconst/+Constant/+QCD/Boltzmann.m b/physconst/+Constant/+QCD/Boltzmann.m
new file mode 100644
index 0000000..88fa71a
--- /dev/null
+++ b/physconst/+Constant/+QCD/Boltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=Boltzmann()
+
+v=1;
+u=9.1e-7;
+
+end
diff --git a/physconst/+Constant/+QCD/LightSpeed.m b/physconst/+Constant/+QCD/LightSpeed.m
new file mode 100644
index 0000000..02df4a3
--- /dev/null
+++ b/physconst/+Constant/+QCD/LightSpeed.m
@@ -0,0 +1,6 @@
+function [v,u]=LightSpeed()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+QCD/Planck.m b/physconst/+Constant/+QCD/Planck.m
new file mode 100644
index 0000000..685f195
--- /dev/null
+++ b/physconst/+Constant/+QCD/Planck.m
@@ -0,0 +1,6 @@
+function [v,u]=Planck()
+
+v=2*pi;
+u=7.5e-8;
+
+end
diff --git a/physconst/+Constant/+QCD/ProtonMass.m b/physconst/+Constant/+QCD/ProtonMass.m
new file mode 100644
index 0000000..cb50ef4
--- /dev/null
+++ b/physconst/+Constant/+QCD/ProtonMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ProtonMass()
+
+v=1;
+u=4.4e-8;
+
+end
diff --git a/physconst/+Constant/+QCD/ReducedPlanck.m b/physconst/+Constant/+QCD/ReducedPlanck.m
new file mode 100644
index 0000000..6792f15
--- /dev/null
+++ b/physconst/+Constant/+QCD/ReducedPlanck.m
@@ -0,0 +1,6 @@
+function [v,u]=ReducedPlanck()
+
+v=1;
+u=1.2e-8;
+
+end
diff --git a/physconst/+Constant/+RA/Boltzmann.m b/physconst/+Constant/+RA/Boltzmann.m
new file mode 100644
index 0000000..88fa71a
--- /dev/null
+++ b/physconst/+Constant/+RA/Boltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=Boltzmann()
+
+v=1;
+u=9.1e-7;
+
+end
diff --git a/physconst/+Constant/+RA/ConductanceQuantum.m b/physconst/+Constant/+RA/ConductanceQuantum.m
new file mode 100644
index 0000000..ee407c4
--- /dev/null
+++ b/physconst/+Constant/+RA/ConductanceQuantum.m
@@ -0,0 +1,6 @@
+function [v,u]=ConductanceQuantum()
+
+v=2/pi;
+u=2e-10;
+
+end
diff --git a/physconst/+Constant/+RA/Coulomb.m b/physconst/+Constant/+RA/Coulomb.m
new file mode 100644
index 0000000..983cead
--- /dev/null
+++ b/physconst/+Constant/+RA/Coulomb.m
@@ -0,0 +1,6 @@
+function [v,u]=Coulomb()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+RA/ElectronMass.m b/physconst/+Constant/+RA/ElectronMass.m
new file mode 100644
index 0000000..81ffd2f
--- /dev/null
+++ b/physconst/+Constant/+RA/ElectronMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ElectronMass()
+
+v=0.5;
+u=2.2e-8;
+
+end
diff --git a/physconst/+Constant/+RA/ElementaryCharge.m b/physconst/+Constant/+RA/ElementaryCharge.m
new file mode 100644
index 0000000..cf86ca9
--- /dev/null
+++ b/physconst/+Constant/+RA/ElementaryCharge.m
@@ -0,0 +1,6 @@
+function [v,u]=ElementaryCharge()
+
+v=sqrt(2);
+u=3.1e-8;
+
+end
diff --git a/physconst/+Constant/+RA/FineStructure.m b/physconst/+Constant/+RA/FineStructure.m
new file mode 100644
index 0000000..a75b283
--- /dev/null
+++ b/physconst/+Constant/+RA/FineStructure.m
@@ -0,0 +1,6 @@
+function [v,u]=FineStructure()
+
+v=7.2973525698e-3;
+u=2.4e-12;
+
+end
diff --git a/physconst/+Constant/+RA/FluxQuantum.m b/physconst/+Constant/+RA/FluxQuantum.m
new file mode 100644
index 0000000..83a291d
--- /dev/null
+++ b/physconst/+Constant/+RA/FluxQuantum.m
@@ -0,0 +1,6 @@
+function [v,u]=FluxQuantum()
+
+v=pi/sqrt(2);
+u=4.9e-8;
+
+end
diff --git a/physconst/+Constant/+RA/LightSpeed.m b/physconst/+Constant/+RA/LightSpeed.m
new file mode 100644
index 0000000..889fa97
--- /dev/null
+++ b/physconst/+Constant/+RA/LightSpeed.m
@@ -0,0 +1,6 @@
+function [v,u]=LightSpeed()
+
+v=274.07199814;
+u=8.8e-8;
+
+end
diff --git a/physconst/+Constant/+RA/Planck.m b/physconst/+Constant/+RA/Planck.m
new file mode 100644
index 0000000..685f195
--- /dev/null
+++ b/physconst/+Constant/+RA/Planck.m
@@ -0,0 +1,6 @@
+function [v,u]=Planck()
+
+v=2*pi;
+u=7.5e-8;
+
+end
diff --git a/physconst/+Constant/+RA/ProtonMass.m b/physconst/+Constant/+RA/ProtonMass.m
new file mode 100644
index 0000000..e922fc9
--- /dev/null
+++ b/physconst/+Constant/+RA/ProtonMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ProtonMass()
+
+v=918.076335;
+u=8.1e-5;
+
+end
diff --git a/physconst/+Constant/+RA/ReducedPlanck.m b/physconst/+Constant/+RA/ReducedPlanck.m
new file mode 100644
index 0000000..6792f15
--- /dev/null
+++ b/physconst/+Constant/+RA/ReducedPlanck.m
@@ -0,0 +1,6 @@
+function [v,u]=ReducedPlanck()
+
+v=1;
+u=1.2e-8;
+
+end
diff --git a/physconst/+Constant/+S/Boltzmann.m b/physconst/+Constant/+S/Boltzmann.m
new file mode 100644
index 0000000..88fa71a
--- /dev/null
+++ b/physconst/+Constant/+S/Boltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=Boltzmann()
+
+v=1;
+u=9.1e-7;
+
+end
diff --git a/physconst/+Constant/+S/Coulomb.m b/physconst/+Constant/+S/Coulomb.m
new file mode 100644
index 0000000..983cead
--- /dev/null
+++ b/physconst/+Constant/+S/Coulomb.m
@@ -0,0 +1,6 @@
+function [v,u]=Coulomb()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+S/ElementaryCharge.m b/physconst/+Constant/+S/ElementaryCharge.m
new file mode 100644
index 0000000..432a486
--- /dev/null
+++ b/physconst/+Constant/+S/ElementaryCharge.m
@@ -0,0 +1,6 @@
+function [v,u]=ElementaryCharge()
+
+v=1;
+u=2.2e-8;
+
+end
diff --git a/physconst/+Constant/+S/Gravitational.m b/physconst/+Constant/+S/Gravitational.m
new file mode 100644
index 0000000..2e47506
--- /dev/null
+++ b/physconst/+Constant/+S/Gravitational.m
@@ -0,0 +1,6 @@
+function [v,u]=Gravitational()
+
+v=1;
+u=4.7e-5;
+
+end
diff --git a/physconst/+Constant/+S/LightSpeed.m b/physconst/+Constant/+S/LightSpeed.m
new file mode 100644
index 0000000..02df4a3
--- /dev/null
+++ b/physconst/+Constant/+S/LightSpeed.m
@@ -0,0 +1,6 @@
+function [v,u]=LightSpeed()
+
+v=1;
+u=0;
+
+end
diff --git a/physconst/+Constant/+SI/AtomicMass.m b/physconst/+Constant/+SI/AtomicMass.m
new file mode 100644
index 0000000..9f19c37
--- /dev/null
+++ b/physconst/+Constant/+SI/AtomicMass.m
@@ -0,0 +1,6 @@
+function [v,u]=AtomicMass()
+
+v=1.660538921e-27;
+u=7.3e-35;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/+SI/AvogadroNumber.m b/physconst/+Constant/+SI/AvogadroNumber.m
new file mode 100644
index 0000000..c9915b9
--- /dev/null
+++ b/physconst/+Constant/+SI/AvogadroNumber.m
@@ -0,0 +1,6 @@
+function [v,u]=AvogadroNumber()
+
+v=6.02214129e23;
+u=2.7e16;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/+SI/BohrMagneton.m b/physconst/+Constant/+SI/BohrMagneton.m
new file mode 100644
index 0000000..686f9b3
--- /dev/null
+++ b/physconst/+Constant/+SI/BohrMagneton.m
@@ -0,0 +1,6 @@
+function [v,u]=BohrMagneton()
+
+v=9.27400968e-24;
+u=2e-31;
+
+end
diff --git a/physconst/+Constant/+SI/BohrRadius.m b/physconst/+Constant/+SI/BohrRadius.m
new file mode 100644
index 0000000..ab6a309
--- /dev/null
+++ b/physconst/+Constant/+SI/BohrRadius.m
@@ -0,0 +1,6 @@
+function [v,u]=BohrRadius()
+
+v=5.2917721092e-11;
+u=1.7e-20;
+
+end
diff --git a/physconst/+Constant/+SI/Boltzmann.m b/physconst/+Constant/+SI/Boltzmann.m
new file mode 100644
index 0000000..f46e725
--- /dev/null
+++ b/physconst/+Constant/+SI/Boltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=Boltzmann()
+
+v=1.3806488e-23;
+u=1.3e-29;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/+SI/ConductanceQuantum.m b/physconst/+Constant/+SI/ConductanceQuantum.m
new file mode 100644
index 0000000..85a1855
--- /dev/null
+++ b/physconst/+Constant/+SI/ConductanceQuantum.m
@@ -0,0 +1,6 @@
+function [v,u]=ConductanceQuantum()
+
+v=7.7480917346e-5;
+u=2.5e-14;
+
+end
diff --git a/physconst/+Constant/+SI/Coulomb.m b/physconst/+Constant/+SI/Coulomb.m
new file mode 100644
index 0000000..4891db4
--- /dev/null
+++ b/physconst/+Constant/+SI/Coulomb.m
@@ -0,0 +1,6 @@
+function [v,u]=Coulomb()
+
+v=8.987551787;
+u=0;
+
+end
diff --git a/physconst/+Constant/+SI/ElectronMass.m b/physconst/+Constant/+SI/ElectronMass.m
new file mode 100644
index 0000000..1390e46
--- /dev/null
+++ b/physconst/+Constant/+SI/ElectronMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ElectronMass()
+
+v=9.10938291e-31;
+u=4e-38;
+
+end
diff --git a/physconst/+Constant/+SI/ElementaryCharge.m b/physconst/+Constant/+SI/ElementaryCharge.m
new file mode 100644
index 0000000..c6bb067
--- /dev/null
+++ b/physconst/+Constant/+SI/ElementaryCharge.m
@@ -0,0 +1,6 @@
+function [v,u]=ElementaryCharge()
+
+v=1.602176565e-19;
+u=3.5e-27;
+
+end
diff --git a/physconst/+Constant/+SI/Faraday.m b/physconst/+Constant/+SI/Faraday.m
new file mode 100644
index 0000000..4375378
--- /dev/null
+++ b/physconst/+Constant/+SI/Faraday.m
@@ -0,0 +1,6 @@
+function [v,u]=Faraday()
+
+v=96485.3365;
+u=2.1e-3;
+
+end
diff --git a/physconst/+Constant/+SI/FermiCoupling.m b/physconst/+Constant/+SI/FermiCoupling.m
new file mode 100644
index 0000000..54d23ad
--- /dev/null
+++ b/physconst/+Constant/+SI/FermiCoupling.m
@@ -0,0 +1,6 @@
+function [v,u]=FermiCoupling()
+
+v=1.166364e-5;
+u=5e-11;
+
+end
diff --git a/physconst/+Constant/+SI/FineStructure.m b/physconst/+Constant/+SI/FineStructure.m
new file mode 100644
index 0000000..a75b283
--- /dev/null
+++ b/physconst/+Constant/+SI/FineStructure.m
@@ -0,0 +1,6 @@
+function [v,u]=FineStructure()
+
+v=7.2973525698e-3;
+u=2.4e-12;
+
+end
diff --git a/physconst/+Constant/+SI/FluxQuantum.m b/physconst/+Constant/+SI/FluxQuantum.m
new file mode 100644
index 0000000..a23671d
--- /dev/null
+++ b/physconst/+Constant/+SI/FluxQuantum.m
@@ -0,0 +1,6 @@
+function [v,u]=FluxQuantum()
+
+v=2.067833758e-15;
+u=4.6e-23;
+
+end
diff --git a/physconst/+Constant/+SI/Gravitational.m b/physconst/+Constant/+SI/Gravitational.m
new file mode 100644
index 0000000..9853df2
--- /dev/null
+++ b/physconst/+Constant/+SI/Gravitational.m
@@ -0,0 +1,6 @@
+function [v,u]=Gravitational()
+
+v=6.67408e-11;
+u=3.1e-15;
+
+end
diff --git a/physconst/+Constant/+SI/HartreeEnergy.m b/physconst/+Constant/+SI/HartreeEnergy.m
new file mode 100644
index 0000000..5b0e999
--- /dev/null
+++ b/physconst/+Constant/+SI/HartreeEnergy.m
@@ -0,0 +1,6 @@
+function [v,u]=HartreeEnergy()
+
+v=4.35974434e-18;
+u=1.9e-25;
+
+end
diff --git a/physconst/+Constant/+SI/LightSpeed.m b/physconst/+Constant/+SI/LightSpeed.m
new file mode 100644
index 0000000..ebaf077
--- /dev/null
+++ b/physconst/+Constant/+SI/LightSpeed.m
@@ -0,0 +1,6 @@
+function [v,u]=LightSpeed()
+
+v=299792458;
+u=0;
+
+end
diff --git a/physconst/+Constant/+SI/Loschmidt.m b/physconst/+Constant/+SI/Loschmidt.m
new file mode 100644
index 0000000..bb1eb39
--- /dev/null
+++ b/physconst/+Constant/+SI/Loschmidt.m
@@ -0,0 +1,6 @@
+function [v,u]=Loschmidt()
+
+v=2.6867805e25;
+u=2.4e19;
+
+end
diff --git a/physconst/+Constant/+SI/NuclearMagneton.m b/physconst/+Constant/+SI/NuclearMagneton.m
new file mode 100644
index 0000000..dee8ef5
--- /dev/null
+++ b/physconst/+Constant/+SI/NuclearMagneton.m
@@ -0,0 +1,6 @@
+function [v,u]=NuclearMagneton()
+
+v=5.05078353e-27;
+u=1.1e-34;
+
+end
diff --git a/physconst/+Constant/+SI/Planck.m b/physconst/+Constant/+SI/Planck.m
new file mode 100644
index 0000000..2e28656
--- /dev/null
+++ b/physconst/+Constant/+SI/Planck.m
@@ -0,0 +1,6 @@
+function [v,u]=Planck()
+
+v=6.626070040e-34;
+u=8.1e-42;
+
+end
diff --git a/physconst/+Constant/+SI/ProtonMass.m b/physconst/+Constant/+SI/ProtonMass.m
new file mode 100644
index 0000000..57163cd
--- /dev/null
+++ b/physconst/+Constant/+SI/ProtonMass.m
@@ -0,0 +1,6 @@
+function [v,u]=ProtonMass()
+
+v=1.672621777e-27;
+u=7.4e-35;
+
+end
diff --git a/physconst/+Constant/+SI/ReducedPlanck.m b/physconst/+Constant/+SI/ReducedPlanck.m
new file mode 100644
index 0000000..74e46cd
--- /dev/null
+++ b/physconst/+Constant/+SI/ReducedPlanck.m
@@ -0,0 +1,6 @@
+function [v,u]=ReducedPlanck()
+
+v=1.054571800e-34;
+u=1.3e-42;
+
+end
diff --git a/physconst/+Constant/+SI/Rydberg.m b/physconst/+Constant/+SI/Rydberg.m
new file mode 100644
index 0000000..36af7a0
--- /dev/null
+++ b/physconst/+Constant/+SI/Rydberg.m
@@ -0,0 +1,6 @@
+function [v,u]=Rydberg()
+
+v=10973731.568539;
+u=5.5e-5;
+
+end
diff --git a/physconst/+Constant/+SI/StefanBoltzmann.m b/physconst/+Constant/+SI/StefanBoltzmann.m
new file mode 100644
index 0000000..7c6891a
--- /dev/null
+++ b/physconst/+Constant/+SI/StefanBoltzmann.m
@@ -0,0 +1,6 @@
+function [v,u]=StefanBoltzmann()
+
+v=5.670373e-8;
+u=2.1e-13;
+
+end
diff --git a/physconst/+Constant/+SI/VacuumImpedance.m b/physconst/+Constant/+SI/VacuumImpedance.m
new file mode 100644
index 0000000..002322e
--- /dev/null
+++ b/physconst/+Constant/+SI/VacuumImpedance.m
@@ -0,0 +1,6 @@
+function [v,u]=VacuumImpedance()
+
+v=376.730313461;
+u=0;
+
+end
diff --git a/physconst/+Constant/+SI/VacuumPermeability.m b/physconst/+Constant/+SI/VacuumPermeability.m
new file mode 100644
index 0000000..d287b25
--- /dev/null
+++ b/physconst/+Constant/+SI/VacuumPermeability.m
@@ -0,0 +1,6 @@
+function [v,u]=VacuumPermeability()
+
+v=1.256637061e-6;
+u=0;
+
+end
diff --git a/physconst/+Constant/+SI/VacuumPermittivity.m b/physconst/+Constant/+SI/VacuumPermittivity.m
new file mode 100644
index 0000000..0ea213e
--- /dev/null
+++ b/physconst/+Constant/+SI/VacuumPermittivity.m
@@ -0,0 +1,6 @@
+function [v,u]=VacuumPermittivity()
+
+v=8.854187817e-12;
+u=0;
+
+end
diff --git a/physconst/+Constant/+SI/VonKlitzing.m b/physconst/+Constant/+SI/VonKlitzing.m
new file mode 100644
index 0000000..27b62a7
--- /dev/null
+++ b/physconst/+Constant/+SI/VonKlitzing.m
@@ -0,0 +1,6 @@
+function [v,u]=VonKlitzing()
+
+v=25812.8074434;
+u=8.4e-6;
+
+end
diff --git a/physconst/+Constant/+SI/WienDisplacement.m b/physconst/+Constant/+SI/WienDisplacement.m
new file mode 100644
index 0000000..58b8e7c
--- /dev/null
+++ b/physconst/+Constant/+SI/WienDisplacement.m
@@ -0,0 +1,6 @@
+function [v,u]=WienDisplacement()
+
+v=2.8977721e-3;
+u=2.6e-9;
+
+end
diff --git a/physconst/+Constant/AtomicMass.m b/physconst/+Constant/AtomicMass.m
new file mode 100644
index 0000000..5ee6f6b
--- /dev/null
+++ b/physconst/+Constant/AtomicMass.m
@@ -0,0 +1,5 @@
+function [v,u]=AtomicMass()
+
+[v,u]=Constant.(UnitSystem()).AtomicMass;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/AvogadroNumber.m b/physconst/+Constant/AvogadroNumber.m
new file mode 100644
index 0000000..0298ad3
--- /dev/null
+++ b/physconst/+Constant/AvogadroNumber.m
@@ -0,0 +1,5 @@
+function [v,u]=AvogadroNumber()
+
+[v,u]=Constant.(UnitSystem()).AvogadroNumber;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/BohrMagneton.m b/physconst/+Constant/BohrMagneton.m
new file mode 100644
index 0000000..5b2d0a1
--- /dev/null
+++ b/physconst/+Constant/BohrMagneton.m
@@ -0,0 +1,5 @@
+function [v,u]=BohrMagneton()
+
+[v,u]=Constant.(UnitSystem()).BohrMagneton;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/BohrRadius.m b/physconst/+Constant/BohrRadius.m
new file mode 100644
index 0000000..d58f519
--- /dev/null
+++ b/physconst/+Constant/BohrRadius.m
@@ -0,0 +1,5 @@
+function [v,u]=BohrRadius()
+
+[v,u]=Constant.(UnitSystem()).BohrRadius;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Boltzmann.m b/physconst/+Constant/Boltzmann.m
new file mode 100644
index 0000000..d626381
--- /dev/null
+++ b/physconst/+Constant/Boltzmann.m
@@ -0,0 +1,5 @@
+function [v,u]=Boltzmann()
+
+[v,u]=Constant.(UnitSystem()).Boltzmann;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/ConductanceQuantum.m b/physconst/+Constant/ConductanceQuantum.m
new file mode 100644
index 0000000..40e45d2
--- /dev/null
+++ b/physconst/+Constant/ConductanceQuantum.m
@@ -0,0 +1,5 @@
+function [v,u]=ConductanceQuantum()
+
+[v,u]=Constant.(UnitSystem()).ConductanceQuantum;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Coulomb.m b/physconst/+Constant/Coulomb.m
new file mode 100644
index 0000000..843552c
--- /dev/null
+++ b/physconst/+Constant/Coulomb.m
@@ -0,0 +1,5 @@
+function [v,u]=Coulomb()
+
+[v,u]=Constant.(UnitSystem()).Coulomb;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/ElectronMass.m b/physconst/+Constant/ElectronMass.m
new file mode 100644
index 0000000..28bd8b0
--- /dev/null
+++ b/physconst/+Constant/ElectronMass.m
@@ -0,0 +1,5 @@
+function [v,u]=ElectronMass()
+
+[v,u]=Constant.(UnitSystem()).ElectronMass;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/ElementaryCharge.m b/physconst/+Constant/ElementaryCharge.m
new file mode 100644
index 0000000..32e4016
--- /dev/null
+++ b/physconst/+Constant/ElementaryCharge.m
@@ -0,0 +1,5 @@
+function [v,u]=ElementaryCharge()
+
+[v,u]=Constant.(UnitSystem()).ElementaryCharge;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Faraday.m b/physconst/+Constant/Faraday.m
new file mode 100644
index 0000000..c8bd741
--- /dev/null
+++ b/physconst/+Constant/Faraday.m
@@ -0,0 +1,5 @@
+function [v,u]=Faraday()
+
+[v,u]=Constant.(UnitSystem()).Faraday;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/FermiCoupling.m b/physconst/+Constant/FermiCoupling.m
new file mode 100644
index 0000000..12c74aa
--- /dev/null
+++ b/physconst/+Constant/FermiCoupling.m
@@ -0,0 +1,5 @@
+function [v,u]=FermiCoupling()
+
+[v,u]=Constant.(UnitSystem()).FermiCoupling;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/FineStructure.m b/physconst/+Constant/FineStructure.m
new file mode 100644
index 0000000..cbd2ac3
--- /dev/null
+++ b/physconst/+Constant/FineStructure.m
@@ -0,0 +1,5 @@
+function [v,u]=FineStructure()
+
+[v,u]=Constant.(UnitSystem()).FineStructure;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/FluxQuantum.m b/physconst/+Constant/FluxQuantum.m
new file mode 100644
index 0000000..eb1dbeb
--- /dev/null
+++ b/physconst/+Constant/FluxQuantum.m
@@ -0,0 +1,5 @@
+function [v,u]=FluxQuantum()
+
+[v,u]=Constant.(UnitSystem()).FluxQuantum;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Gravitational.m b/physconst/+Constant/Gravitational.m
new file mode 100644
index 0000000..e3f56c2
--- /dev/null
+++ b/physconst/+Constant/Gravitational.m
@@ -0,0 +1,5 @@
+function [v,u]=Gravitational()
+
+[v,u]=Constant.(UnitSystem()).Gravitational;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/HartreeEnergy.m b/physconst/+Constant/HartreeEnergy.m
new file mode 100644
index 0000000..3cb8b49
--- /dev/null
+++ b/physconst/+Constant/HartreeEnergy.m
@@ -0,0 +1,5 @@
+function [v,u]=HartreeEnergy()
+
+[v,u]=Constant.(UnitSystem()).HartreeEnergy;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/LightSpeed.m b/physconst/+Constant/LightSpeed.m
new file mode 100644
index 0000000..b6e2357
--- /dev/null
+++ b/physconst/+Constant/LightSpeed.m
@@ -0,0 +1,5 @@
+function [v,u]=LightSpeed()
+
+[v,u]=Constant.(UnitSystem()).LightSpeed;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Loschmidt.m b/physconst/+Constant/Loschmidt.m
new file mode 100644
index 0000000..f09764e
--- /dev/null
+++ b/physconst/+Constant/Loschmidt.m
@@ -0,0 +1,5 @@
+function [v,u]=Loschmidt()
+
+[v,u]=Constant.(UnitSystem()).LoschmidtConstant;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/NuclearMagneton.m b/physconst/+Constant/NuclearMagneton.m
new file mode 100644
index 0000000..4b50b21
--- /dev/null
+++ b/physconst/+Constant/NuclearMagneton.m
@@ -0,0 +1,5 @@
+function [v,u]=NuclearMagneton()
+
+[v,u]=Constant.(UnitSystem()).NuclearMagneton;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Planck.m b/physconst/+Constant/Planck.m
new file mode 100644
index 0000000..9021d6a
--- /dev/null
+++ b/physconst/+Constant/Planck.m
@@ -0,0 +1,5 @@
+function [v,u]=Planck()
+
+[v,u]=Constant.(UnitSystem()).Planck;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/ProtonMass.m b/physconst/+Constant/ProtonMass.m
new file mode 100644
index 0000000..a979470
--- /dev/null
+++ b/physconst/+Constant/ProtonMass.m
@@ -0,0 +1,5 @@
+function [v,u]=ProtonMass()
+
+[v,u]=Constant.(UnitSystem()).ProtonMass;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/ReducedPlanck.m b/physconst/+Constant/ReducedPlanck.m
new file mode 100644
index 0000000..fa4974a
--- /dev/null
+++ b/physconst/+Constant/ReducedPlanck.m
@@ -0,0 +1,5 @@
+function [v,u]=ReducedPlanck()
+
+[v,u]=Constant.(UnitSystem()).ReducedPlanck;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/Rydberg.m b/physconst/+Constant/Rydberg.m
new file mode 100644
index 0000000..5f0f8d7
--- /dev/null
+++ b/physconst/+Constant/Rydberg.m
@@ -0,0 +1,5 @@
+function [v,u]=Rydberg()
+
+[v,u]=Constant.(UnitSystem()).Rydberg;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/StefanBoltzmann.m b/physconst/+Constant/StefanBoltzmann.m
new file mode 100644
index 0000000..607cb4f
--- /dev/null
+++ b/physconst/+Constant/StefanBoltzmann.m
@@ -0,0 +1,5 @@
+function [v,u]=StefanBoltzmann()
+
+[v,u]=Constant.(UnitSystem()).StefanBoltzmann;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/VacuumImpedance.m b/physconst/+Constant/VacuumImpedance.m
new file mode 100644
index 0000000..e76b379
--- /dev/null
+++ b/physconst/+Constant/VacuumImpedance.m
@@ -0,0 +1,5 @@
+function [v,u]=VacuumImpedance()
+
+[v,u]=Constant.(UnitSystem()).VacuumImpedance;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/VacuumPermeability.m b/physconst/+Constant/VacuumPermeability.m
new file mode 100644
index 0000000..ae42938
--- /dev/null
+++ b/physconst/+Constant/VacuumPermeability.m
@@ -0,0 +1,5 @@
+function [v,u]=VacuumPermeability()
+
+[v,u]=Constant.(UnitSystem()).VacuumPermeability;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/VacuumPermittivity.m b/physconst/+Constant/VacuumPermittivity.m
new file mode 100644
index 0000000..cf82b5c
--- /dev/null
+++ b/physconst/+Constant/VacuumPermittivity.m
@@ -0,0 +1,5 @@
+function [v,u]=VacuumPermittivity()
+
+[v,u]=Constant.(UnitSystem()).VacuumPermittivity;
+
+end
\ No newline at end of file
diff --git a/physconst/+Constant/WienDisplacement.m b/physconst/+Constant/WienDisplacement.m
new file mode 100644
index 0000000..edc7622
--- /dev/null
+++ b/physconst/+Constant/WienDisplacement.m
@@ -0,0 +1,5 @@
+function [v,u]=WienDisplacement()
+
+[v,u]=Constant.(UnitSystem()).WienDisplacement;
+
+end
\ No newline at end of file
diff --git a/physconst/README.md b/physconst/README.md
new file mode 100644
index 0000000..410e5e5
--- /dev/null
+++ b/physconst/README.md
@@ -0,0 +1,63 @@
+# PhysConst
+A collection of fundamental physical constants in various unit systems for MATLAB.
+
+## Licensing
+This software is licensed under the GNU General Public License (version 3).
+
+## Tested On
+- MATLAB R2013b - R2018a
+
+## Syntax
+`v=Constant.NameOfConstant` returns the value of the named constant in the default unit system. See the [list of physical constant names](https://github.com/heriantolim/PhysConst#list-of-physical-constant-names). For examples:
+- `v=Constant.ElementaryCharge` returns the value of the elementary charge, *e*;
+- `v=Constant.LightSpeed` returns the value of the speed of light, *c*.
+
+### Default Unit System
+By default, the values are returned in the [SI](https://en.wikipedia.org/wiki/International_System_of_Units) units. To set a different unit system as default, do one of the following:
+- Make a variable assignment: `UnitSystem='UnitSystemCode';`, where `UnitSystemCode` refers to the unit system to be used. See the [list of unit system codes](https://github.com/heriantolim/PhysConst#list-of-unit-system-codes).
+- Define an anonymous function: `UnitSystem=@()'UnitSystemCode';`.
+- Create a function named `UnitSystem` in the current working directory or elsewhere. The function takes no input arguments and returns a string specifying the unit system code.
+
+`v=Constant.UnitSystemCode.NameOfConstant` returns the value of the named constant in the specified unit system, regardless of the default unit system. For examples:
+- `v=Constant.SI.ElementaryCharge` returns the value of the elementary charge, *e*, in the SI units (=1.602176565);
+- `v=Constant.P.LightSpeed` returns the value of the speed of light, *c*, in the Planck units (=1).
+
+### Uncertainties
+`[v,u]=Constant.NameOfConstant` and `[v,u]=Constant.UnitSystemCode.NameOfConstant` additionally return the uncertainty of the named constant.
+
+## List of Unit System Codes
+- `HA`: Hartree Atomic Units
+- `P`: Planck Units
+- `QCD`: Quantum Chromodynamics Units
+- `RA`: Rydberg Atomic Units
+- `S`: Stoney Units
+- `SI`: International System of Units
+
+## List of Physical Constant Names
+- `AtomicMass`: atomic mass unit, *m*<sub>u</sub>
+- `AvogadroNumber`: Avogadro's number, *N*<sub>A</sub>
+- `BohrMagneton`: Bohr magneton, *μ*<sub>B</sub>
+- `BohrRadius`: Bohr radius, *a*<sub>0</sub>
+- `Boltzmann`: Boltzmann constant, *k*<sub>B</sub>
+- `ConductanceQuantum`: conductance quantum, *G*<sub>0</sub>
+- `Coulomb`: Coulomb constant, *k*<sub>e</sub>
+- `ElectronMass`: electron mass, *m*<sub>e</sub>
+- `ElementaryCharge`: elementary charge, *e*
+- `Faraday`: Faraday constant, *F*
+- `FermiCoupling`: Fermi coupling constant, *G*<sub>F</sub>/(*ħc*)<sup>3</sup>
+- `FineStructure`: fine-structure constant, *α*
+- `FluxQuantum`: magnetic flux quantum, Φ<sub>0</sub>
+- `Gravitational`: Newtonian constant of gravitation, *G*
+- `HartreeEnergy`: Hartree energy, *E*<sub>h</sub>
+- `LightSpeed`: speed of light, *c*
+- `Loschmidt`: Loschmidt's number, *n*<sub>0</sub>
+- `NuclearMagneton`: nuclear magneton, *μ*<sub>N</sub>
+- `Planck`: Planck constant, *h*
+- `ProtonMass`: proton mass, *m*<sub>p</sub>
+- `ReducedPlanck`: reduced Planck constant, *ħ*
+- `Rydberg`: Rydberg constant, *R*<sub>∞</sub>
+- `StefanBoltzmann`: Steffan-Boltzmann constant, *σ*
+- `VacuumImpedeance`: vacuum impedance, *Z*<sub>0</sub>
+- `VacuumPermeability`: vacuum permeability, *μ*<sub>0</sub>
+- `VacuumPermittivity`: vacuum permittivity, *ε*<sub>0</sub>
+- `WienDisplacement`: Wien displacement law constant, *b*
diff --git a/physconst/UnitSystem.m b/physconst/UnitSystem.m
new file mode 100644
index 0000000..24dbe92
--- /dev/null
+++ b/physconst/UnitSystem.m
@@ -0,0 +1,5 @@
+function s=UnitSystem()
+
+s='SI';
+
+end
\ No newline at end of file