VerilogA issues -- transition function

My code is below. I was hesitant to post it because it's such a niche application but the model is for a ferroelectric capacitor. The model has weird current spikes at times however and also fails to converge at times, which is what I'm trying to fix. The spikes in my opinion occur when polarization_next and polarization_prev have the same value on and off for some time span (those are the curves I posted initially) and I was wondering if it would be possible to impose a condition such that the polarization_next and polarization_prev variable would not be evaluated at such frequent time intervals. Perhaps this would get rid of the plateaus but I don't know whether this is possible in VerilogA. I was experimenting with the timer function for this purpose.

I'm wondering if the following line of code is wrong which is what I'm investigating right now. Unfortunately it's difficult to work with 1D arrays in verilogA with something that should be a 2D array....

domain_states[j] = 1;

I will get back to you with the above line if I find that it's wrong.

There are some papers on the topic like "Preisach model for the simulation of ferroelectric capacitors" by Bartic and " Preisach Polarization - Electric Field Hysteresis model with iteration method" by Hideyuki Ikeda.

Basically the model calculates the polarization (also the name of the variable) for an applied voltage V(p,n). The current is proportional to the time derivative of the polarization. I drew a sketch that hopefully helps but I'm not sure how much it will help without reading any of the papers. The polarization is calculated by multiplying each element in the domain_states array by each corresponding element in the preisach_density array and summing the results.

Thank you for any help.

##################################################################### code starts here ###########################################

/*

This file implements the extraction of csv data into a verilog variable and then implements a Preisach implemenatation of a FE capacitor SPICE
model based on the Preisach data.

*/

`include "constants.vams"
`include "disciplines.vams"

module test_func(p,n);

    inout p,n;
    electrical p,n;

    //files for reading/writing
    integer mcd, fp;

    real data, v_present, v_step;

    real preisach_density[0:1024];// all elements should add up to 1
    real domain_states[0:1024];// all elements should be -1 or +1 (depending on polarity)
    real polarization;//polarization (between -Pr and +Pr)
    
    real polarization_prev, polarization_next;
    real time_prev, time_next;
    real v_prev, v_next;

    integer row_offset, column_offset;

    real current_calc;

    integer i, j, k;//used in loops

    integer rows, columns;//describes area scanned from applied voltage
    integer rows_total, columns_total;//gives number of rows and columns in the grid ("rows" and "columns" should not exceed these values)
    
    analog begin    

        @(initial_step) begin

            //read data from csv file and transfer it to verilog-A array variable
            //data stored in csv file should be preisach distribution
            //csv file should have the same number of elements as arrays "preisach_density" and "domain_states"
            mcd=$fopen("whole_preisach_dist_5_1_2017.csv","r");

            while (!$feof(mcd)) begin

                $fscanf(mcd,"%f", preisach_density[i]);
                domain_states[i] = 0;

                i = i+1;
                  end

        v_step = 0.1; // resolution

        rows_total = 3.2/v_step;//the maximum anticipated applied input voltage swing divided by the resolution of the grid cells
        columns_total = 3.2/v_step;//the maximum anticipated applied input voltage swing divided by the resolution of the grid cells

        v_next = V(p,n);

        time_prev = 0;
        time_next = 0;
        
        // centers the indiciess to the beginning of the fourth quadrant of the Preisach distribution (i.e. half of row_total or columns total plus 1)
        row_offset = 17;
        column_offset = 17;
        
        end



        if (v_next < v_prev) begin // negative voltage derivative

            rows = row_offset + V(p,n)/(-v_step);
            
            begin
            if(rows > rows_total)
                rows = rows_total;
            else if(rows < 0)
                rows = 0;
            end

            //for all columns
            for(j = 0; j < columns_total; j=j+1) begin

                //for all rows with a smaller value than "rows"
                for(k=0; k<rows; k=k+1) begin

                    domain_states[k*columns_total + j] = -1;
                end
            end
        end

        if (v_next > v_prev) begin // positive voltage derivative

            columns = column_offset + V(p,n)/(v_step);

            begin
            if(columns > columns_total)
                columns = columns_total;
            else if(columns < 0)
                columns = 0;
            end

            //value compared to j should equal number of columns in array
            for (j=0; j<columns * rows_total ; j=j+1) begin

                domain_states[j] = 1;

            end
        end

            //calculate the present polarization across the device
            polarization = 0;
            for (i=0; i<columns_total*rows_total; i = i+1) begin

                polarization = polarization + preisach_density[i]*domain_states[i];
            end
            
            //used for calculating derivatives
            polarization_prev = polarization_next;
            polarization_next = polarization;
            v_prev = v_next;
            v_next = V(p,n);
            time_prev = time_next;
            time_next = $abstime;

            begin
                if (time_next == 0)

                    current_calc = 0;

                else
                    current_calc = (polarization_next - polarization_prev)/(time_next - time_prev)*0.000000001;//0.000000001 is arbitrary

            end
        
        I(p,n) <+ current_calc;
        
    end

endmodule

My code is below. I was hesitant to post it because it's such a niche application but the model is for a ferroelectric capacitor. The model has weird current spikes at times however and also fails to converge at times, which is what I'm trying to fix. The spikes in my opinion occur when polarization_next and polarization_prev have the same value on and off for some time span (those are the curves I posted initially) and I was wondering if it would be possible to impose a condition such that the polarization_next and polarization_prev variable would not be evaluated at such frequent time intervals. Perhaps this would get rid of the plateaus but I don't know whether this is possible in VerilogA. I was experimenting with the timer function for this purpose.

I'm wondering if the following line of code is wrong which is what I'm investigating right now. Unfortunately it's difficult to work with 1D arrays in verilogA with something that should be a 2D array....

domain_states[j] = 1;

I will get back to you with the above line if I find that it's wrong.

There are some papers on the topic like "Preisach model for the simulation of ferroelectric capacitors" by Bartic and " Preisach Polarization - Electric Field Hysteresis model with iteration method" by Hideyuki Ikeda.

Basically the model calculates the polarization (also the name of the variable) for an applied voltage V(p,n). The current is proportional to the time derivative of the polarization. I drew a sketch that hopefully helps but I'm not sure how much it will help without reading any of the papers. The polarization is calculated by multiplying each element in the domain_states array by each corresponding element in the preisach_density array and summing the results.

Thank you for any help.

##################################################################### code starts here ###########################################

/*

This file implements the extraction of csv data into a verilog variable and then implements a Preisach implemenatation of a FE capacitor SPICE
model based on the Preisach data.

*/

`include "constants.vams"
`include "disciplines.vams"

module test_func(p,n);

    inout p,n;
    electrical p,n;

    //files for reading/writing
    integer mcd, fp;

    real data, v_present, v_step;

    real preisach_density[0:1024];// all elements should add up to 1
    real domain_states[0:1024];// all elements should be -1 or +1 (depending on polarity)
    real polarization;//polarization (between -Pr and +Pr)
    
    real polarization_prev, polarization_next;
    real time_prev, time_next;
    real v_prev, v_next;

    integer row_offset, column_offset;

    real current_calc;

    integer i, j, k;//used in loops

    integer rows, columns;//describes area scanned from applied voltage
    integer rows_total, columns_total;//gives number of rows and columns in the grid ("rows" and "columns" should not exceed these values)
    
    analog begin    

        @(initial_step) begin

            //read data from csv file and transfer it to verilog-A array variable
            //data stored in csv file should be preisach distribution
            //csv file should have the same number of elements as arrays "preisach_density" and "domain_states"
            mcd=$fopen("whole_preisach_dist_5_1_2017.csv","r");

            while (!$feof(mcd)) begin

                $fscanf(mcd,"%f", preisach_density[i]);
                domain_states[i] = 0;

                i = i+1;
                  end

        v_step = 0.1; // resolution

        rows_total = 3.2/v_step;//the maximum anticipated applied input voltage swing divided by the resolution of the grid cells
        columns_total = 3.2/v_step;//the maximum anticipated applied input voltage swing divided by the resolution of the grid cells

        v_next = V(p,n);

        time_prev = 0;
        time_next = 0;
        
        // centers the indiciess to the beginning of the fourth quadrant of the Preisach distribution (i.e. half of row_total or columns total plus 1)
        row_offset = 17;
        column_offset = 17;
        
        end



        if (v_next < v_prev) begin // negative voltage derivative

            rows = row_offset + V(p,n)/(-v_step);
            
            begin
            if(rows > rows_total)
                rows = rows_total;
            else if(rows < 0)
                rows = 0;
            end

            //for all columns
            for(j = 0; j < columns_total; j=j+1) begin

                //for all rows with a smaller value than "rows"
                for(k=0; k<rows; k=k+1) begin

                    domain_states[k*columns_total + j] = -1;
                end
            end
        end

        if (v_next > v_prev) begin // positive voltage derivative

            columns = column_offset + V(p,n)/(v_step);

            begin
            if(columns > columns_total)
                columns = columns_total;
            else if(columns < 0)
                columns = 0;
            end

            //value compared to j should equal number of columns in array
            for (j=0; j<columns * rows_total ; j=j+1) begin

                domain_states[j] = 1;

            end
        end

            //calculate the present polarization across the device
            polarization = 0;
            for (i=0; i<columns_total*rows_total; i = i+1) begin

                polarization = polarization + preisach_density[i]*domain_states[i];
            end
            
            //used for calculating derivatives
            polarization_prev = polarization_next;
            polarization_next = polarization;
            v_prev = v_next;
            v_next = V(p,n);
            time_prev = time_next;
            time_next = $abstime;

            begin
                if (time_next == 0)

                    current_calc = 0;

                else
                    current_calc = (polarization_next - polarization_prev)/(time_next - time_prev)*0.000000001;//0.000000001 is arbitrary

            end
        
        I(p,n) <+ current_calc;
        
    end

endmodule